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Abstract 

A  class  of  multilevel  methods  for  second  order  problems  is  considered  in  the  additive 
Schwarz  framework.  It  is  established  that,  in  the  general  case,  the  condition  number  of 
the  iterative  operator  grows  at  most  linearly  with  the  number  of  levels.  The  bound  is 
independent  of  the  mesh  sizes  and  the  number  of  levels  under  a  regularity  assumption. 
This  is  an  improvement  of  a  result  by  Dryja  and  Widlund  on  a  multilevel  additive  Schwarz 
algorithm,  and  the  theory  given  by  Bramble,  Pasciak  and  Xu  for  the  BPX  algorithm. 

Additive  Schwarz  and  iterative  substructuring  algorithms  for  the  biharmonic  equation 
are  also  considered.  These  are  domain  decomposition  methods  which  have  previously  been 
developed  extensively  for  second  order  elliptic  problems  by  Bramble,  Pasciak  and  Schatz, 
Dryja  and  Widlund  and  others. 

Optimal  convergence  properties  are  established  for  additive  Schwarz  algorithms  for 
the  biharmonic  equation  discretized  by  certain  conforming  finite  elements.  The  number 
of  iterations  for  the  iterative  substructuring  methods  grows  only  as  the  logarithm  of  the 
number  of  degrees  of  freedom  associated  with  a  typical  subregion.  It  is  also  demonstrateed 
that  it  is  possible  to  simplify  the  basic  tdgorithms.  This  leads  to  a  decrease  of  the  cost 
but  not  of  the  rate  of  convergence  of  the  iterative  methods.  In  the  ancilysis,  new  tools  are 
developed  to  deal  with  Hermitian  elements.  Certain  new  inequalities  for  discrete  norms 
for  finite  element  spaces  are  also  used. 
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Chapter  1 

Introduction 

1.1      An  Overview 

Finite  element  discretizations  of  elliptic  problems,  e.g.  the  biharmonic  problem,  often  re- 
sult in  very  large,  sparse  linear  systems.  Domain  decomposition  methods  are  very  powerful 
iterative  methods  for  solving  such  problems.  The  domain  is  decomposed  into  overlapping 
(resp.  non-overlapping)  subdomains;  the  algorithms  are  often  referred  to  as  Schwarz  type 
(resp.  iterative  substructuring)  methods.  In  each  iteration  step,  a  number  of  smaller  sub- 
problems,  which  correspond  to  the  restriction  of  the  original  problem  to  subregions,  or  to 
a  problem  with  a  coarse  mesh,  are  solved.  The  number  of  subproblems  can  be  large  and 
these  methods  are  promising  for  parallel  computation. 

The  Schwarz  alternating  method,  the  first  domain  decomposition  algorithm,  was  pro- 
posed by  H.A.  Schwarz  [42]  in  1869.  It  involves  solving  problems  on  two  subdomains 
sequentially.  In  [34],  P.  L.  Lions  interpreted  the  method  in  a  variational  framework. 
Many  subdomain  cases  were  also  considered.  Because  of  the  sequential  nature  of  the 
method,  it  is  not  ideal  for  parallel  computations.  In  [24],  Dryja  and  Widlund  introduced 
the  additive  Schwarz  method  to  remove  this  sequential  behavior  of  the  Schwarz  alternating 
method  making  it  possible  to  solve  all  the  subproblems  in  parallel,  see  also  Matsokin  and 
Nepomnyaschikh  [37]  and  Nepomnyaschikh  [39].  The  additive  Schwarz  method  can  be 
regarded  as  an  iterative  method  for  an  equivalent  system  with  a  better  condition  number. 
Dryja  and  Widlund  showed  the  optimal  convergence  properties  of  the  algorithm  for  second 
order  problems  discretized  by  linear  elements.    Generalization  to  certain  stationary  and 


parabolic  convection-diffusion  problems  have  been  carried  out  by  Cai  and  Widlund;  cf. 
[14,16,17]. 

Iterative  substructuring  methods  for  second  order  problems  were  developed  by  Bram- 
ble, Pasciak  and  Schatz  [10,11],  Bjorstad  and  Widlund  [7],  Widlund  [47,46]  and  Dryja 
[23,22].  Techniques  for  analyzing  some  iterative  substructuring  methods  using  the  addi- 
tive Schwarz  framework  were  developed  in  Dryja  and  Widlund  [25].  They  demonstrated 
that  certain  iterative  substructuring  methods  can  be  viewed  as  additive  Schwarz  methods 
with  a  set  of  special  overlapping  subdomains.  For  discussions  of  the  relationship  between 
the  two  type  schemes,  see  also  Bjorstad  and  Widlund  [8]  and  Chan  and  Goovaerts  [19]. 

Multilevel  methods,  such  as  multigrid  methods,  are  among  the  most  efficient  methods 
for  linear  equations  arising  from  elliptic  problems;  cf.  Hackbusch  [32],  McCormick  [38] 
and  the  references  therein.  Recently,  with  the  increasing  interest  in  parallel  computation, 
several  new  multilevel  methods  have  been  developed;  cf.  Yserentant  [52],  Bank,  Dupont 
and  Yserentant  [1],  Bramble,  Pasciak  and  Xu  [13],  and  Dryja  and  Widlund  [28].  In  this 
thesis,  we  give  improved  results  for  a  class  of  multilevel  methods  by  showing  that  the 
condition  number  of  the  iteration  operator  grows  at  most  linearly  with  the  number  of 
levels  in  general,  and  is  bounded  by  a  constant  independent  of  the  mesh  sizes  and  the 
number  of  levels  if  the  elliptic  problem  is  fl^^-regular.  This  is  an  improvement  on  Dryja 
and  Widlund's  results  on  a  multilevel  additive  Schwarz  method  as  well  as  Bramble,  Pasciak 
and  Xu's  results  on  the  BPX  algorithm. 

We  also  study  the  additive  Schwarz  and  iterative  substructuring  methods  for  the  bi- 
harmonic  Dirichlet  problem.  We  construct  additive  Schwarz  algorithms  for  the  biharmonic 
equation  using  standard  conforming  finite  element  discretization.  Some  new  tools  are  de- 
veloped for  the  proofs  of  the  optimal  convergences  of  the  algorithms.  By  using  a  weak 
coupling  property  of  the  degrees  of  freedom,  some  simplified  versions  are  also  derived. 
This  leads  to  a  decrease  of  the  cost  per  iteration  but  not  of  the  rate  of  convergence  of  the 
iterative  methods.  In  the  case  of  the  iterative  substructuring  methods,  we  demonstrate 
that  direct  generalizations  of  some  known  iterative  substructuring  methods  result  in  algo- 
rithms with  a  condition  number  which  grows  at  least  like  0((77//i)^).   Better  algorithms 


are  obtained  by  using  certain  vertex  spaces.  Some  parallel  multilevel  algorithms  are  also 
constructed  for  the  biharmonic  problem.  Some  earlier  works  on  the  biharmonic  problem 
can  be  found  in  Glowinski  and  Pironneau  [30],  Bjorstad  [2],  WidJund  [45]  and  Chan,  E 
and  Sun  [18]. 

The  thesis  is  organized  as  follows.  In  the  remainder  of  chapter  1,  we  review  some  basic 
Sobolev  spaces.  We  also  discuss  some  preliminary  material  for  the  biharmonic  problem 
and  some  standard  finite  element  discretizations.  In  chapter  2,  we  review  some  domain 
decomposition  techniques  in  a  form  used  in  this  thesis.  In  chapter  3,  we  discuss  a  class 
of  multilevel  methods,  including  a  multilevel  additive  Schwarz  method  of  Dryja  and  Wid- 
lund  and  the  BPX  algorithm  of  Bramble,  Pasciak  and  Xu.  We  present  our  algorithms  and 
analysis  using  the  additive  Schwarz  framework.  In  chapter  4,  we  construct  some  interpo- 
lation and  quasi-interpolation  operators  and  establish  approximation  properties  of  these 
operators.  We  also  construct  discrete  Sobolev  norms  and  establish  their  equivalence  to  the 
corresponding  continuous  Sobolev  norms.  We  also  prove  a  theorem  regarding  the  norm 
estimates  of  some  standard  interpolation  operators.  These  results  are  used  in  chapters  5 
and  6  to  show  the  optimal  and  almost  optimal  convergence  properties  of  the  algorithms.  In 
chapter  5,  we  construct  additive  Schwarz  schemes  for  the  biharmonic  problem  discretized 
by  some  standard  finite  elements  and  we  establish  the  optimal  convergence  properties  of 
the  algorithms.  By  using  the  weak  coupling  property  of  the  degrees  of  freedom  of  the 
finite  elements,  some  computationally  more  efficient  algorithms  are  also  derived.  We  also 
consider  a  multilevel  method.  Numerical  experiments  with  the  algorithms  are  reported.  In 
chapter  6,  we  discuss  some  generalization  of  the  iterative  substructuring  methods  for  the 
biharmonic  problem.  We  first  demonstrate  that  direct  generalizations  of  some  algorithms 
designed  for  the  second  order  problems  result  in  algorithms  with  very  large  condition 
numbers.  We  then  show  that  better  algorithms  can  be  obtained  by  adding  certain  vertex 
spaces  to  the  space  decomposition.  Finally,  we  present  numericaJ  experiments  to  support 
our  theoretical  conclusions. 


1.2      Some  Sobolev  Spaces 

Let  fi  C  7J''  be  a  bounded  Lipschitz  domain  and  let  L''{Q)  be  the  Banach  space  of  p-th 
power  integrable  functions. 

We  use  the  multi-index  notation  for  partial  derivatives.  Let  o  =  (qi.  •  •  • ,  Od)  be  a  non- 
negative  multi-index  representing  the  order  of  the  partial  derivatives,  |q|  h  qi  -f 1-  qj. 

and  let 

^„4ef     d''     def  d- 


dx"        dx^'-'-dx'^" 
The  Sobolev  space  W^-^iQ)  is  defined  by 

Vr'"'P(fi)  =^  {v\  D^v  6  L^in),  if  \o\  <  m} 
with  a  norm 

\a\<m 

We  also  need  the  semi  norm 

|a|=m 

The  most  useful  space  in  this  thesis  is  the  case  p  =  2,  which  by  convention  is  denoted  by 

To  study  the  homogeneous  Dirichlet  problem,  we  need  the  Sobolev  space 

n^iQ)^^  {v\v  e  ^'"(n),Z)°r|an  =  0,|a|  <  m  -  1}. 
This  is  the  closure  of  the  space  of  Co°(n)  with  respect  to  //'"-norm. 

The  following  inequalities  establish  equivalences  of  certain  norms  in  subspaces  of 

Lemma  1.2.1  (Friedrichs'  inequality)   There  exists  a  positive  constant  C{il),  which 
depends  only  on  the  Lipschitz  constant  of  the  boundary  ofQ,  such  that,  for  all  u  G  HQ{Ct), 

M^ii)  <  C(n)//n|«lw.(n)- 
Here  Hq  is  the  diameter  ofQ. 


Let  {u}n  =  "Wrr  be  the  average  value  of  the  function  u  over  il.  We  have 

Lemma  1.2.2  (Poincar^'s  inequality)   There  exists  a  positive  constant  C{il),   which 
depends  only  on  the  Lipschitz  constant  of  the  boundary  ofil,  such  that,  for  all  u  €  H^(U), 

h  -  {«}nll//'(n)  <  C{Q)Ho\u\ffi(n)- 
Here  Hq  is  the  diameter  ofil. 

Proofs  of  Lemma  1.2.1  and  Lemma  1.2.2  can  be  found  in  [40]. 

For  applications  to  the  biharmonic  problem,  the  following  generalization  of  Poincare's 
inequality  is  useful. 

Lemma  1.2.3  (Generalized  Poincare's  Inequality)  //{Z?'*u}n  =  0,  |qI  <  k  then, 

Mh'IH)  <  Cn^n"^'"'l"l//*+i(n)    fors<k+l,ue  H''-^\il), 

and 

The  second  inequality  also  holds  for  the  W'''^^'^-norm,  with  d/2  replaced  by  d/p. 

Proof.    The  first  inequality  can  be  proved  by  applying  Poincare's  inequality  several 
times  to  D°u.  The  second  follows  from  the  embedding  theorem: 

n^-^\Q)^W''°°{n),     s<k  +  l-d/2; 

the  powers  of  Hq  are  obtained  by  a  scaling  argument.  I 

We  aJso  need  the  following  lemma. 

Lemma  1.2.4  Ifu  €  n'''^^{il),  then  there  exists  a  polynomial  p  €  Vk  such  that 

{Z?>-p)}n  =  0,      |a|<;:. 

Proof.  We  establish  that  for  any  given  da,  \a\  <  k,  there  exists  a  polynomial  p  €  Vk, 
such  that 

I  D''p=da. 


The  lemma  then  follows  by  setting  da  =  Jq  D°u.  Let 

We  define  a  partial  order  for  the  set  of  multi-index  q  =  (01,02,  •  •  -,0^): 

a<0      iffa.  <^.,      l<i<d. 
We  have 

D°p  =     Yl    h^'''^"  =     E    hCa,0T'^-°,  (1.2) 

\P\<k  \0\<k 

where  Ca,p  =  j^^  for  /?  >  a,  and  Co,/?  =  0  otherwise.  We  obtain  a  linear  system 

Y  ^o,0b0^  [  D''p  =  da,  (1.3) 

with  Aq,0  =  Cq^0 ^^x^~°dx.  It  is  easy  to  see  that  Aa,0  =  0  for  /?  ^  q,  and  Aa,a  =  a!|fi|- 
We  can  therefore  find  a  permutation  so  that  the  matrix  Aa,0  is  upper  triangular  with 
nonzero  diagonals.  Therefore  there  exists  an  unique  solution  6^  and  thus  a  unique  p  £  Vk 
satisfying  {D''(u  -  p)}n  =  0,Vue /r*+'(n).  ■ 

Corollary  1.2.5  For  the  polynomial  p  constructed  above,  we  have 

I"  -  p\h-{Q)  <  <^n^n'*"'~*l«l//*+'(n)    ^^^    s<k+l 
Lemma  1.2.6  (Quotient  space  lemma)  Let 

II"IIhvp*-i  =  jgf  ll«-pll//»- 

Pfc^k-l 

Then 

l"l//»  <  ll«ll//vn-r 
For  a  bounded  Lipschtz  region  il,  we  have 

ll«ll//Vp»->  ^  cin)\u\„. 

Proof.  The  lemma  follows  from  Corollary  1.2.5  and  Lemma  1.2.3.  ■ 


Lemma  1.2.7  (Bramble-Hilbert's  lemma)  Let  Q  be  an  open  Lipschitz  region  in  R^. 
Let  k  >  0  be  an  integer  and  p  6  [0,oo].  Let  f  be  a  continuous  linear  form  on  iy*+''''(Q) 
with  the  property  that 

/(p)  =  o,    VpeP*(fi). 

Then,  there  exists  a  constant  C{il)  such  that 

where  ||  •  |lvv*+i,prn)  **  '^^  norm  of  the  dual  space  o/ VF*'*'^'''(fi). 
Proof.  For  all  p  G  /'^(fi),  we  have 

\fiv)\  =  \fiv  +  p)\  <  ||/|IIv*+.,P(n)k  +  Plvy^+..P(n), 

and  thus 

\fiv)\  <  ll/IIIv*+>.p(n)     Jnf  ,  1^  +  Plw*+>.p(n)- 
p€P*(n) 

The  conclusion  follows  by  the  quotient  space  lemma.  I 

1.3     The  Biharmonic  Equation 

Consider  the  biharmonic  Dirichlet  problem  in  a  plane  region 

A^u    =     /         in  n, 

u    =    go        on         dQ,  (1.4) 

^     =     gi        on         dil. 

There  are  also  other  types  of  boundary  conditions  which  are  of  interest.  For  example,  one 
can  impose  u|an  =  9o  and  Au  =  52;  this  corresponds  to  the  mathematical  model  for  a 
simply  supported  plate.  That  problem  can  be  decomposed  into  two  second  order  problems 
and  is  therefore  relatively  easy  to  solve. 

Weak  Formulation.  We  consider  the  variational  form  of  the  problem  with  homogeneous 
boundary  condition:  Find  u  £  ^o(^)  such  that 

a{u,v)  =  f{v),     VveH^iil),  (1.5) 


where  /  is  a  bounded  linear  functional  on  H^i^l)  and  a(u,  v)  is  a  symmetric,  continuous, 
T/Q-elliptic  bilinear  form.  Two  examples  of  such  bilinear  forms  are 

a{u,v)-   /  AuAv  dx,  (1.6) 

and 

a(u,v)=  /  {AuAv+    1-<T    2^-^— ;r— ^ ^-7^^-^-2TT)}<'^'         (1-7 

^n  0x10x20x10x2      dx{dx^      0x2  oxf 

The  first  one  arises  in  Fluid  Dynamics,  and  the  second  provides  a  variational  formulation 
of  the  Clamped  Plate  Problem.  Here  0  <  cr  <  1/2  is  Poisson's  coefficient  of  the  plate.  By 
the  Lax-Milgram  Theorem,  the  boundness  and  ellipticity  imply  existence  and  uniqueness 
of  the  solution;  cf.  Lax  and  Milgram  [33]  and  Ciarlet[21]. 

1.3.1      The  Finite  Element  Formulation 

The  finite  element  formulation  is  obtained  by  replacing  the  infinite  dimensional  space 
V  =  Hq{Q)  with  a  finite  dimensional  subspace  V^  C  V. 

We  triangulate  the  domain  fi  into  non-overlapping  regions  called  elements,  generally 
triangles  or  rectangles.  V''  is  a  space  of  piecewise  polynomials  with  respect  to  the  trian- 
gulation.  The  finite  element  solution  u^  G  V'*  satisfies 

aiuK,<i>h)  =  fi4>h),     V<^,  GV\  (1.8) 

We  note  that  the  finite  element  solution  Uh  is  the  a(-,  •)-projection  of  the  true  solution  u 
onto  the  finite  element  space  V'*. 

The  norm  derived  from  the  bilinear  form  a(-,-)  defines  a  semi-norm  in  the  Sobolev 
space  H^{^).  It  is  a  norm  of  the  space  HQ{il),  and  therefore  a  norm  in  its  subspace  V^. 

Let  {</>,}  be  the  nodal  basis  for  V''.  Then  Uh  can  be  represented  as 

I 

Thus  we  obtain  a  linear  system  for  i,  the  degrees  of  freedom  of  Uh, 

Kx  =  b 


with  the  stiffness  matrix  given  by 

and  the  load  vector  by 

The  stiffness  matrix  K  is  symmetric,  positive  definite.  After  a  proper  scaling,  its  con- 
dition number  k{K)  =  0{h~*).  Since  the  system  is  usually  very  large,  and  the  condition 
number  of  K  is  also  quite  large,  solving  the  system  can  be  very  expensive.  Many  precon- 
ditioners  have  been  designed  for  K.  Among  them,  the  additive  Schwarz  methods  studied 
in  this  thesis  seem  to  be  particularly  successful  and  promising. 

1.3.2     Some  Conforming  Elements 

For  the  biharmonic  equation,  the  finite  elements  are  all  relatively  complicated.  In  this 
thesis,  we  restrict  ourselves  to  some  standard  conforming  elements.  In  particular,  we 
consider  the  Argyris  triangle  V^,  the  Bell  triangle  V^  and  the  bicubic  element  Vq.  These 
elements  are  complicated  but  among  the  simplest  conforming  elements  for  the  biharmonic 
equation. 

The  Argyris  element  consists  of  continuously  differentiable  functions,  the  restrictions 
of  which  to  any  triangle  are  in  V5.  The  degrees  of  freedom  for  this  element  in  a  triangle 
with  vertices  ai,i  =  1,2,3,  are  given  by 

d"  d 

where  bi  is  the  midpoint  of  edge  ajait,  and  n,  is  the  outward  normal  direction  ofajok.  The 
number  of  the  degrees  of  freedom  for  each  triangle  is  21. 

It  is  easy  to  see  that,  in  generad,  the  normal  derivatives  of  an  Argyris  element  is  a 
polynomial  of  degree  4.  Let  Vb  be  the  subspace  of  V5  formed  by  those  polynomials  of 
V5  whose  normal  derivatives  along  each  side  of  a  triangle  are  polynomials  of  degree  3  in 
t,  the  abscissa  along  an  axis  containing  the  side.  We  note  that  V4  C  Vb  C  P5.  The  Bell 
element  consists  of  C^  functions  whose  restrictions  to  a  triangle  are  in  Vb-  The  degrees 


of  freedom  for  the  Bell  element  are  given  by 


{p{a.),£p(a.),|^p(a.),£^p(a.),^;^p(«.),^p(a.)}- 


a.\ 
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Figure  1.1:  The  Argyris  and  Bell  triangles 

If  the  domain  is  build  from  rectangles,  we  can  also  use  the  bicubic  element,  known  as 
the  Bogner-Fox-Schmit  rectangle  in  the  engineering  literature.  It  is  space  of  C'  functions 
whose  restrictions  to  a  rectangle  are  in  Q3  =  span {I'y-',  i,j  <  3}.  The  degrees  of  freedom 
of  the  bicubic  element  are  given  by 

We  note  that  nonconforming  elements  such  as  Morley's  triangle,  Adini's  rectangle, 
etc.,  are  also  widely  used  in  engineering  computation.  However,  they  will  not  be  discussed 
here. 

The  approximation  properties  of  these  elements  are  well  understood.  We  have 

||«-«^||2.n<C'/i^|ti|6.n, 
||«-UB||2.n<C/i^l«|5.n. 
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Figure  1.2:  The  bicubic  element 

II"-  "Qslb.n  <  Ch'^\u\4^Q, 
For  more  details  on  the  convergence  and  error  estimates;  see  Ciarlet  [20,21]. 

1.3.3     Properties  of  the  Basis  Functions 

Only  the  following  two  properties  of  the  basis  functions  will  be  used. 

•  The  uniformity  of  the  basis  functions. 

•  The  invariance  of  linear  functions  under  interpolation. 

We  state  the  properties  for  the  Bell  element,  the  same  results  hold  for  the  Argyris  and 
bicubic  elements.  Let  us  denote  the  nodal  basis  functions  of  the  Bell  element  by  <pf,  \a\  < 
2.  The  basis  functions  satisfy 


dx0 


ar     >.      /  1,     ifi  =  i 
'^^  (^^^  =  \  0,      otherw 


=  t  and  /?  =  q; 
erwise. 


Lemma  1.3.1    The  basis  functions  <f>f  of  the  Bell  element  are  uniform  to  order  2  in  the 
sense  of  Strang  [44],  »•«• 

|<^rk'.~  <C'/i''*l-',5<2,  (1.9) 

and  uniform  to  order  2  in  the  following  weak  sense, 

\4>i\H'  <C/i'+l"l-%5<2. 
The  same  conclusions  hold  for  the  Argyris  element  and  the  bicubic  element. 
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(1.10) 


If  we  write  down  the  basis  functions  explicitly,  it  is  straight  forward  to  check  that  they 
are  uniform  to  order  2. 

We  remark  that  for  the  additive  Schwarz  methods,  we  only  need  ( 1.10).  For  some  basis 
functions,  it  is  easy  to  check  that  (1.10)  holds,  but  difficult  to  check  whether  (1.9)  holds 
or  not. 

Another  property  of  the  basis  functions  is  the  invariance  of  linear  function  under 
interpolation.  Let  ^^/^  be  the  standard  interpolation  operator  to  V^.  For  p  linear  in  a 
triangle,  Uyhp  =  p.  This  implies  that,  restricted  to  a  triangle, 

p  G  span{0f;|a|  <  1}, 

which  gives  us  the  following  relationships  for  the  basis  functions  of  the  Bell  element: 

1    =    E.</'.(x) 

y   =   E.y.'^.(^)  +  <^f(a^) 

Similar  equations  can  be  worked  out  for  the  Argyris  elements  and  the  bicubic  elements. 

The  above  equations  show  that  in  each  triangle,  a  linear  function  can  be  constructed 
from  only  the  basis  functions  associated  with  nodaJ  values  and  the  first  derivatives,  since 
the  basis  functions  associated  with  second  derivatives  do  not  appear  in  the  above  equation. 
This  fact  is  important  in  later  chapters. 
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Chapter  2 

Domain  Decomposition 
Techniques 

2.1      The  Conjugate  Gradient  Method 

Stationary  schemes  like  the  Jacobi,  Gauss-Seidel,  SOR,  Block  Jacobi  and  Block  Gauss- 
Seidel  methods  can  be  viewed  as  subspace  correction  methods.  The  approximate  solution  is 
updated  by  solving  problems  associated  with  a  set  of  subspaces  of  the  solution  space.  These 
methods  can  be  accelerated  by  using  the  conjugate  gradient  method  or  the  Chebyshev 
semi-iterative  method;  we  concentrate  on  the  conjugate  gradient  method. 

The  conjugate  gradient(CG)  method  is  of  fundamental  importance  for  the  domain  de- 
composition methods.  All  algorithms  studied  in  this  thesis  are  variations  of  preconditioned 
conjugate  gradient  methods(PCG).  Since  CG  and  PCG  are  well  known  algorithms,  we  only 
give  a  brief  description.  For  more  complete  treatment,  see  Golub  and  Van  Loan  [31].  We 
note  that  in  the  conjugate  gradient  algorithm,  we  only  need  to  be  able  to  apply  the  matrix 
j4  to  a  given  vector;  an  explicit  representation  of  the  matrix  is  not  needed.  This  is  a  very 
important  for  domain  decomposition  algorithms. 

Conjugate  Gradient  Algorithm 

Set  A:  =  0;  lo  =  0;  To  =  b. 
while  |rfc|  >  f|ro| 
Jt  =  Jt-|-l 
if  A;=  1 

Pi  =  ro 
else 

13k  =  ('•fc_i,rjt_i)/(rfc_2,rjt-2) 
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Pk  -  rh-\  +  0kPk-i 
end 

OJt  =  irk,rk)/(Pk,Apk) 
Xk  =  Xk  +  ajtp* 
r*:  =  rt  -  QkAp'' 
end 

The  reduction  in  the  energy  norm  of  the  error,  after  n  steps  of  conjugate  gradient  iteration, 
is  given  by,  see  e.g.  [31], 

lk-^n||<2(^^r|k-Zoll. 

Here  k  =  k{A)  is  the  condition  number  of  A  given  by 

_   "m&x(A) 
^min(A) 

The  conjugate  gradient  method  is  therefore  often  an  effective  iterative  algorithm  to  solve 
symmetric,  positive  definite  systems 

Ax  =  b. 

The  algorithm  converges  to  a  good  approximate  solution  in  relatively  few  iterations  for  a 
well  conditioned  matrix  A.  When  A  is  not  well  conditioned,  which  is  generally  the  case 
for  discretizations  of  elliptic  problems,  we  can  introduce  a  preconditioner  B  and  solve  the 
preconditioned  linear  system 

B-^Ax  =  B-^b. 

by  using  the  conjugate  gradient  method  with  the  i4-inner  product.  Domain  decomposition 
is  an  effective  way  to  construct  such  preconditioners.  Other  preconditioners  are  also  used 
in  practices,  e.g.  the  diagonal  scaling  and  incomplete  Cholesky  preconditioners. 

2.2     Multiplicative  Schwarz  Schemes 

The  Schwarz  alternating  method  is  the  first  domain  decomposition  method,  proposed 
by  H.  A.  Schwarz  in  1869,  [42].  Schwarz  established  convergence  using  the  maximum 
principle.  In  the  1930's,  Sobolev  extended  the  result  to  the  partial  differential  equations 
of  linear  elasticity  [43].  New  theoretical  tools  for  the  analysis  of  multiplicative  Schwarz 
schemes  have  been  developed  more  recently.  The  continuous  problem  has  been  studied  by 
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Lions  [34];  see  also  Nepomnyaschikh  [39].  For  the  finite  element  case,  the  two  subspace 
case  has  been  discussed  in  Bjorstad  [3],  Mandel  and  McCormick  [35]  and  Bjorstad  and 
Mandel  [4].  Results  for  the  case  of  more  then  two  subspaces  have  been  found  by  Widlund 
[48]  and  Mathew  [36]  and  most  successfully  by  Bramble,  Pasciak,  Wang,  and  Xu  [12].  We 
review  this  theory  and  the  results  which  can  be  used  in  analyzing  multiplicative  variants 
of  the  algorithms. 

Let  us  consider  the  domain,  shown  in  Figure  2.1,  with  fi  =  fii  U  fij  on  which  we  wish 

to  solve 

r  -Au  =  /  in  fi, 

1       u  =  0  on  dQ. 


Figure  2.1:  Schwarz  alternating  scheme 


The  Schwarz  alternating  procedure  approximates  the  solution  iteratively  by  solving 
problems  on  the  each  subdomain  ft,. 


and 


Let 


i 
{ 


_^yn+l/2  ^  y 
„n+l/2  ^  „n 


in  fii, 
on  dQi, 

in  ^2, 
on  0^2- 


a(u,  v)=  /  Vu  •  Vr       and       /(f)  =  /  fv. 
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Then  the  algorithm  can  be  written  as 

a(u"+'/2  -  „",</.)  =  f{<t>)  -  a(u",<^),  «"+i/2  -u-^e  i/J(0,),V<A  6  H^Q,), 

and 

In  each  half-step,  the  correction  is  thus  the  projection  of  the  error  onto  the  subspace 
H^iili)  or  ^^(fia).  The  projection  P,  :  H^iQ)  -»  H^{n,),  is  defined  by 

The  error  propagation  operator  for  a  complete  step  of  the  alternating  Schwarz  method  is 
simply 

{I  -  P2)il  -  Pi). 

The  alternating  Schwarz  method  can  be  considered  as  a  special  example  of  an  abstract 
multiplicative  Schwarz  scheme.  Consider  the  general  variational  problem:  Find  u  E  V 
such  that 

a{u,4>)  =  fi<f>),     V<^€V.  (2.1) 

Let  Vi  be  subspaces  of  V  so  that  V  =  Vi  + \-  Vjv-  Associated  with  each  subspace  is  the 

corresponding  orthogonaJ  projection  operator  P,  :  V  — ►  V,.  The  abstract  multiplicative 
Schwarz  scheme  proceeds  by  sequentially  projecting  the  error  onto  the  subspaces  Vj.  The 
error  propagation  operator  for  a  complete  step  of  the  alternating  Schwarz  method  is  given 

by 

{l-p^)...{i-p,). 

» 

To  include  the  case  of  approximate  solver  for  the  subproblems,  it  is  sometimes  necessary 
to  consider  a  more  general  form 

(/_r^)...(/_T,), 

where  T,  are  symmetric,  positive  semi-definite,  linear  operators  with 

||r.lU<u;<2. 
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Usually,  the  T,  are  either  the  projections  P,  or  approximations  thereof.  Define  Ei  by 

iE;.  =  (/-T.)(/-T,_i).-.(/-Ti). 

The  effectiveness  of  the  particular  multiplicative  algorithm  is  determined  by  the  reduction 
in  the  error  for  one  complete  iteration  step,  i.e.  the  best  7  in  the  inequality 


\ElMa  <  l\v\ 


Estimates  for  7  are  given  in  Bramble,  Pasciak,  Wang  and  Xu  [12]  and  Xu  [51]. 
Let 

and  let  XauniT)  be  the  minimum  eigenvalue  of  T.    A  main  result  of  Bramble,  Pasciak, 
Wang,  and  Xu  is  given  in 
Theorem  2.2.1 

\Env\1  <  ^\v\l 

leith 

,      (2-u;)A^n(r)  {2-u)\r^{T) 

7  -        max^^^  ^  a;2Ar(A^  _  l)/2)'  2(1  +  u^N{N  -  l)/2y' 

To  obtain  stronger  results,  we  use  a  coarse  subspace  Vq  =  V^  and  make  assumptions 
on  the  subspaces  Vj.  Let  €  be  the  matrix  such  that 

\<Vi,  vj)\  <  eija{vi,  VifMvj.  Vj)'f\  Vt;^  G  K,  \/vj  €  V,,  t,;  =  1,  •  •  •  iV.      (2.2) 

The  inequalities  (2.2)  are  strengthened  Cauchy-Schwarz  inequalitie^s.  Let  p(S)  be  the 
spectral  radius  of  S.  We  note  that  as  the  overlap  between  the  spaces  increases  Amin(^) 
will  generally  improve  while  the  bound  in  (2.2)  will  degrade,  i.e.  the  spectral  radius  p{S) 
will  increase. 

Using  the  strengthened  Cauchy  inequality  (2.2),  the  above  result  can  be  improved;  cf. 
Xu  [51]. 
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Theorem  2.2.2  Assume  that  (2.2)  hold.   Then 

\Env\1  <  l\v\l 

with 

^  (2  -  u>)An^„(T) 

''  2(l+u;V(^)2)' 

An  important  parameter  in  the  above  two  theorems  is  Xmn{T).  Estimates  for  XtmniT)  was 
obtained  by  Dryja  and  Widlund  [24,25].  An  important  special  case  of  this  genera]  result 
is  when  the  T,  are  projections  F,.  In  this  case,  \mn{T)  can  be  estimated  by  Cq^  in  the 
following  assumption:  For  all  r  G  V^,  there  exists  a  decomposition  f  =  X)  ^>  ^'^^  ^<  ^  ^ 

such  that 

N 

^a(t;,,t;,)  <  Cla{v,v). 
1=1 

For  the  proof  of  this  result,  see  Lemma  2.3.2  in  the  following  section.  In  addition,  we  note, 
that  since  the  Pi  are  orthogonal  projections,  u;  is  one. 

We  remark  that  it  is  also  possible  to  treat  one  or  several  of  the  subspaces  separately; 
cf.  Bramble,  Pasciak,  Wang  and  Xu  [12],  Dryja  [27]  and  Widlund  [49]. 

2.3      Additive  Schwarz  Schemes 

One  possible  problem  with  the  multiplicative  Schwarz  methods  is  the  sequential  nature 
of  the  fractional  steps  of  each  iteration.  For  more  than  two  subspaces,  these  schemes  are 
also  inherently  nonsymmetric.  Therefore,  symmetric  variants,  which  can  almost  double 
the  number  of  fractional  steps,  have  to  be  used  if  we  wish  to  use  the  standard  conjugate 
gradient  method  to  accelerate  the  convergence.  The  additive  Schwarz  methods  were  de- 
signed by  Dryja  and  Widlund  to  remove  the  inherent  sequential  behavior  of  the  fractional 

* 

steps  and  to  preserve  symmetry.  Independent  work  on  additive  Schwarz  methods  can 
also  be  found  in  Matsokin  and  Nepomnyaschikh  [37]  and  Nepomnyaschikh  [39].  Optimal 
convergence  properties  of  certain  algorithms  were  established  for  second  order  self-adjoint 
elliptic  problems,  see  [24,25].  Generalizations  to  nonsymmetric  and  indefinite  cases  have 
been  made  by  Cai  and  Widlund;  cf.[14,15,16,17].  We  will  show  in  chapter  5  that  certain 
additive  Schwarz  methods  for  the  biharmonic  equation  are  also  optimal. 
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Additive  Schwarz  methods  can  be  viewed  as  iterative  methods  for  the  solution  of  an 
auxiliary  linear  problem  that  has  the  same  solution  as  the  original  finite  element  problem. 
We  consider  the  finite  dimensional  variational  problem:  Find  u  £  \'  such  that 

a{u,4>)  =  fief,),     W<t>eV.     '  (2.3) 

Suppose  V  can  be  written  as  the  sum  of  iV  +  1  subspaces 

Let  P,  :  V  — *  Vi,i  -  0,  •  •  • ,  iV,  be  the  orthogonal  projections,  with  respect  to  a(-,  •),  defined 
by 

a{P,u,<f>)  =  a{u,4>),     ^4>eV,.  (2.4) 

Let  P  =  ^Pi-  The  additive  Schwarz  method  for  solving  equation  (2.3  i  is  introduced  in 
terms  of  an  auxiliary  problem:  Find  u  E  V  by  solving 

Pu  =  J2  Ptu  =  g.  (2.5) 

i 

The  right  hand  side  g  has  to  be  chosen  such  that  the  auxiliary  equation  has  the  same 
solution  as  equation  (2.3).  For  this  to  hold,  the  right-hand  side  must  be  equal  to  g  = 
^i=o9ii  where  j,  =  P.u/,  can  be  computed  by  solving 

a{g„v)  =  a{PiUh,v)  =  a{uh,v)  =  f{v),\/v  e  V;. 

Thus  PiU  can  be  found  without  knowing  the  solution  of  (2.3).  Once  the  right-hand  side 
g  is  known,  we  can  use  an  iterative  method,  e.g.  the  conjugate  gradient  method,  to  solve 
equation  (2.5).  To  define  an  additive  Schwarz  method,  it  is  sufficient  to  define  the  space 
decomposition.  The  fundamental  issue  i»  to  estimated  the  condition  number  of  P. 

The  reason  for  going  from  problem  (2.3)  to  problem  (2.5)  is  that,  by  a  suitable  choice 
of  the  subspaces  Vi,  we  can  turn  a  large  ill-conditioned  system  into  a  very  well  conditioned 
problem  at  the  expense  of  solving  many  small  independent  linear  systems. 

The  following  lemmas,  which  establish  the  relations  between  the  decomposition  of  u 
and  the  spectrum  of  P,  allow  us  to  develop  bounds  on  the  condition  number  of  P. 
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Lemma  2.3.1   Let  V  be  a  Hilbert  space,  Vi  be  subspaces  ofV  and  V  =  J2^i-  ^^'  A',  be 
the  projections  from  V  to  V,  and  P  =  Yli  Pv,-   Then  P  is  invertible  and 

aiP-^  u,  u)  =    min    ^  a{ui,  u,) 

with  u,  6  V,  and  the  minimum  is  achieved  at  u,  =  P,P~^u. 

Proof.   We  use  the  properties  of  the  projections  and  Cauchy-Schwarz's  inequality  to 
obtain, 

=     j:a{P,P-^u,u,) 

=  a(p-'u,uy/HLMiy^'. 

Thus, 

a{P~^u,u)  <  ^a{u„Ui). 

Equality  holds  if  and  only  if  u,  =  P,P~^u.  I 

The  following  lemma  is  a  direct  consequence  of  Lemma  2.3.1. 

Lemma  2.3.2  Let  V  be  a  Hilbert  space,  V,  be  subspaces  ofV  and  V  =  YL^i-  ^^^  Pv,  ^ 
the  projections  from  V  to  Vi  and  P  =  Yli  Pv,  ■   Then 

-^minC^)  = '^max(/'     )  =  max —  =  max    min    ■ — 

"        a[u,u)  "    E"'="      «(">") 

and 

^mt^iP)  =  >^Tmn{P     )  =  "iin —  =  mm    nun ^— . 

ti        a{u.u)  u    ^u,=u      a(u,u) 

Remark  2.3.1  If  we  can  find  a  constant  C\  such  that  there  exists  a  decomposition  of 
u  =  J2t  "i  satisfying 

^a(u,,  u.)  <  Cia(u,u),     Vu  G  V, 

then,  it  follows  from  Lemma  2.3.2  that  Anii„(P)  >  C^^.  This  result,  known  as  Lions' 
lemma,  is  very  important  in  estimating  the  minimum  eigenvalue  of  P\  cf.  Dryja  and 
Widlund  [24]. 
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If  we  can  find  a  constant  C2  such  that  Vu  G  F  and  for  any  decomposition  of  u  =  Yli  ^ii 
we  have 

^a(u,,u,)  >  C2a(u,u) 

then,  from  the  Lemma  2.3.2,  we  know  that  XmaxiP)  <  Cf  ^. 

For  domain  decomposition  algorithms,  the  following  lemma,  gives  a  convenient  way  of 
estimating  an  upper  bound  of  A(P). 

Lemma  2.3.3  Consider  the  undirected  graph  with  a  node  for  each  subspace  V,,  and  an 
edge  between  node  i  and  node  j  if  and  only  «/ V,  n  Vj  ^  0.  Let  Nc  be  the  number  of  colors 
needed  to  color  the  graph  so  that  no  two  nodes  connected  by  an  edge  have  the  same  color. 
Then 

Azn»x(/')  <  Nc. 

Proof.  All  the  subspaces  of  a  particular  color  are  disjoint;  hence  their  corresponding 
projection  operators  are  mutually  orthogonal.  Therefore  the  sum  of  the  projection  opera- 
tors of  a  particular  color  is  itself  a  projection  operator.  P  then  is  the  sum  of  Nc  projection 
operators,  each  of  norm  one.  I 

An  alternative  way,  as  in  the  multiplicative  case,  is  to  use  a  strengthened  Cauchy- 
Schwarz  inequality.  The  proof  of  the  following  lemma  is  quite  elementary. 

Lemma  2.3.4  Assume  equation  (2.2)  holds,  then 
This  is  especially  useful  for  the  multilevel  schemes. 

» 

Remark  2.3.2  As  in  the  multiplicative  case,  we  can  also  use  approximate  solvers  T,  to 
replace  P,. 

Remark  2.3.3  In  the  case  of  two  subspaces,  the  relation  between  the  additive  and  mul- 
tiplicative versions  for  the  same  pair  of  subspaces  is  well  understood;  cf.  [3], [4],  and  [35]. 
However  little  can  be  said  for  the  case  of  more  than  two  subspaces.  The  same  two  pa- 
rameters p{€)  and  AnunC^)  provide  bounds  for  the  convergence  rates  of  both  variants. 
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However,  there  is  no  general  theory  which  explicitly  relates  the  actual  convergence  rates 
of  the  two  variants.  The  architectures  of  a  particular  parallel  machine,  etc.  seems  more 
likely  to  determine  which  variant  to  prefer  rather  than  a  strict  mathematical  analysis. 
We  note  that  the  multiplicative  variant  leads  to  a  nonsymmetric  operator.  Therefore  to 
accelerate  it  with  the  conjugate  gradient  method,  we  must  introduce  additional  fractional 
steps  to  make  it  symmetric.  We  could  also  accelerate  the  nonsymmetric  problem  directly 
using  GMRES,  a  conjugate  gradient  type  method;  cf.  Eisenstat,  Elman  and  Schultz  [29] 
and  Saad  and  Schultz  [41]. 
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Chapter  3 

Second  Order  Problems 

3.1      Second  Order  Equations 

We  consider  second  order  uniformly  elliptic  equations  of  the  following  type 

.^(a,(.)^)u     =     /     in       fi,  ^^^^ 

u     =     g     on    dil, 

with  {flij}  symmetric,  uniformly  positive  definite,  bounded  and  piecewise  smooth  on  a 
Lipschitz  region  CI  C  R'^,  d  =  2  or  3. 

It  is  convenient  to  write  (3.1)  in  a  weak  form:  Find  u  6  ^o(")  such  that 

aiu,v)^f{v),^veH^iQ),  (3.2) 

with  the  bilinear  form  a{u,  v)  and  the  functional  f{v)  given  by 

du      dv 


dx. 


The  bilinear  form  a(u,  v)  is  bounded  and  uniformly  elliptic  (coercive),  i.e.    there  exist 
constants  Ci  and  C2  such  that 


a(u,t;)<Ci|«|;,^(n)l^l//^(n)       and       a(u,u)  >  CjIu]^^ 

Existence  and  uniqueness  follows  from  the  Lax-Milgram's  Theorem. 
To  simplify  the  presentation,  we  work  with  the  model  problem 

f  -Au(a:)     =     /(i)     in       fi, 
I  u{x)    =        0       on     5fi, 
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(«) 


(3.3) 


We  consider  a  finite  element  discretization  of  equation  (3.2).  We  triangulate  the  do- 
main n  into  non-overlapping  elements  r,  (triangles  or  rectangles)  and  denote  the  triangu- 
lation  by  T^  =  {r,}.  We  assume  the  triangulation  is  shape  regular.  The  finite  element 
space  V^  is  a  subspace  of  ^o(^)  consisting  of  continuous  functions  whose  restrictions  to 
each  element  r,  are  polynomials.  Here,  we  will  restrict  ourselves  to  the  simplest  finite 
elements,  in  particular,  we  will  use  the  linear  element  (for  R^  or  R^),  the  bi-linear  element 
(for  7?^)  or  the  tri-linear  element  (for  R^).  We  will  use  the  nodal  basis  functions  {<t>,{x)} 
as  the  basis  for  V''.  The  nodal  basis  functions  satisfy 

,^.(1,)  =  6,„     Vx,  6  A\fi). 

To  obtain  the  finite  element  problem,  we  replace  the  Sobolev  space  HqI^)  by  the  finite 
dimensional  subspace  V'':  Find  Uh  6  V*  such  that 

aiuH,VH)  =  f{vH),     Vrh6F\  (3.4) 

Using  the  nodal  basis  {((>,},  u(x)  can  be  represented  as  u{x)  =  ^,x,4>,(x).  Thus,  we 
obtain  a  linear  system 

Kx  =  6,  (3.5) 

where  K  =  {a{(f>i,(pj)]  is  the  stiffness  matrix  and  b,  =  f(4>i)-  ^  is  a  symmetric,  positive 

definite,  sparse  matrix.  It  is  known  that  its  condition  number  k{A)  =  0{h~^).  Thus  A"  is 

■ 
ill  conditioned  for  large  problems. 

3.2     Additive  Schwarz  Methods 
3.2.1      The  Algorithm 

From  chapter  2,  we  know  that  in  order  to  define  an  additive  Schwarz  scheme,  it  is  sufficient 
to  define  the  space  decomposition;  different  space  decompositions  give  different  algorithms. 
Following  Dryja  and  Widlund  [24,25],  we  define  two  levels  of  triangulations.  We  start 
with  a  coarse  triangulation  T^  =  {fi,}.  Each  fi,  is  then  further  divided  into  smaller 
elements  to  get  a  fine  triangulation  T^  =  {r,}.  Let  ^,=diameter  of  fi,,  H  =  max,{F,} 
/ii=diameter  of  r,,  and  h  =  majCi{/i,}. 
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To  get  an  overlapping  decomposition  of  fi,  we  extend  each  fi,  to  a  larger  region  fi,, 
such  that  cHi  <  dist{d^i,dCl,}  <  CHi,  and  dCli  align  with  the  boundaries  of  element  r,. 
We  cut  off  the  part  of  (),  that  is  outside  of  Q. 
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Figure  3.1:  Two  levels  of  triangulation 

It  is  easy  to  see  that  {fi,}  forms  a  finite  covering  of  domain  Q.  Due  to  the  generous 
overlap  of  {(li),  we  have  a  partition  of  unity  {^,}  satisfying 

N 

^e,  =  l  with  0i  e  W^'°^{Cl,),0<  e,  <  l,  and  I^.Ih-i.oo  <  C/H,. 
1=1 

Let  V''  and  V^  be  the  linear  elements  associated  with  the  triangulations  T^  and  T^, 
respectively.  Let  Vo  =  V"  and  Vi  =  V'^{Cli)  =  V'^  n  n^{^i).  We  obtain  a  decomposition 

of  the  finite  element  space  V'' 

N 

t=0 

Let  Pv,  :  V^  — >  Vi,  be  the  ^^-projection  defined  by 

a{Pv,u,v)  =  a{u,v),     Vu  G  K. 
The  matrix  form  of  Py, ,  after  a  permutation  is 


^'■  =  (o''''   SJ'^^^M-. 
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where  K  is  the  stiffness  matrix  associated  with  the  domain  Q  and  A",  is  the  stiffness  matrix 
associated  with  the  subdomain  fi,.  The  additive  Schwarz  operator,  P  is  given  by 

N 

1=0 

The  additive  Schwarz  algorithm  for  equation  (3.4)  is  to  solve  an  auxilary  equation  which 
is  equivalent  to  equation  (3.4). 

Algorithm  3.2.1  Find  u/,  E  V^  by  solving  itemtively  the  equation 

Puh  =  9h,  (3.6) 

with  gh  =  J2i  9x  ^^^  '^^  S'  given  by  the  solutions  for  the  following  problems 

a{gi,4>h)  =  a{P,u,<t>h)  =  fiMMh  e  v;.  (3.7) 

We  use  the  conjugate  gradient  method  to  solve  Pu^  =  g^.  We  only  need  to  compute  Pv 
for  a  given  r  G  F''  in  each  iteration,  thus,  the  explicit  representation  of  P  is  not  needed. 
The  computation  is  carried  out  in  the  following  steps. 

•  Compute  the  right  hand  side  g^  by  first  solving  equations  (3.7)  and  then  setting 
9  =  12  9i- 

•  In  each  iteration,  compute  Pvh  for  given  Vh  by  first  solving 

a{PiVh,4>)  =  aivh,<f>),     '^4>eVi, 

and  then  setting  Pv^^  =  ^Z  PiV^ 

Remark  3.2.1  We  consider  a  special  choice  of  subregions.  Let  J),  =  supp{<;!!),},  and  thus 
V,  =  span{</),}.  Then  the  corresponding  additive  Schwarz  method,  without  a  coarse  space, 
corresponds  to  applying  the  CG  algorithm  to  D~^Kx  =  D~^b,  where  D  —  diag{K).  This 
is  the  Jacobi  conjugate  gradient  method. 

3.2.2     Condition  Number  Estimate 

Theorem  S.2.1  (Dryja  and  Widlund)  For  the  operator  P  defined  above,  there  exists 
a  constant  C\  such  that 

Cia{v,v)<  a{Pv,v)<  (Nc+  l)a{v,v),     Vv  6  V^. 
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Thus  k{P)  <  (Nc  +  1)/Ci,  and  the  rate  of  convergence  of  the  additive  Schwarz  algorithm 
is  independent  of  the  mesh  size  and  the  number  of  subdomains. 

Proof.  We  first  establish  the  upper  bound.  It  is  trivial  to  see  that 

a(Pu,u)  <  (A^  +  l)a(u,u). 

However  we  want  to  get  an  upper  bound  independent  of  the  number  of  substructures  as 
well  as  the  mesh  size  h.  For  u  G  V^, 

a{PiU,u)  =  a{PiU,Piu)  =  a^  (P,u,P,u)  <  a^  (k,u). 

Summing  over  j,  using  the  finite  covering  property  of  {fl,},  we  get 

N 
y^  a{PiU,  u)  <  Nca{u,  u). 

:=1 

We  also  note  that  a{Pou,  u)  <  a{u,  u).  Thus  the  upper  bound  follows  with  C2  =  {Nc  +  1). 

To  prove  the  lower  bound,  for  any  u  G  V^'',  one  needs  to  find  a  good  partition  {«,}  of 
u  as  in  Lions'  Lemma.  Let  Q^  :  ^o(^)  ~*  ^^  be  the  L2  projection.  It  is  known  that  Q^ 
has  the  following  property 

\\^-Q"^\\mo)    <    C^KI/ZMn)  (3-8) 

\\Q"^\\hhu)    <    C\u\H,^o)  (3.9) 

Let  n''  :  ^o(^)  ^  C'(fi)  —*■  V'^i  be  the  standard  interpolation  operator.  Let  uq  =  Ufj  = 
Q^Uh,  Wh  —  Uh  -  Uf{i  and  u,  =  n'^{0iWh).  It  is  easy  to  see  that  Uh  =  J2o  '"■i-  To  see 
that  we  indeed  get  a  good  partition  of  Uh,  we  have  to  estimate  the  norms  of  u,.  It  can  be 
shown  that 

|n''(^.«')l//i(ft.)    <    <^l».k~l«'l//>(ft.)  +  ^l^''^''~l"''i'(n.) 
Using  the  approximation  property  of  Q^ ,  we  have 
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Summing  over  t,  and  using  the  finite  covering  property  of  {n,},  we  obtain 

i 

The  lower  bound  of  P  follows  from  Lions'  Lemma.  I 

Remark  3.2.2  If  we  do  not  use  the  coarse  space,  the  condition  number  of  the  operator 
will  grows  like  n~^.  This  fact  was  pointed  out  and  proved  by  Widlund;  cf.  [47]. 

In  practice,  the  maximum  eigenvalue  of  P  is  between  Nc  and  Nc  +  I,  and  gets  closer 
to  Nc  as  H/h  grows. 

3.3     Multilevel  Methods 

In  two  level  methods,  we  need  to  solve  a  coarse  problem  of  size  0(1/ H"^)  and  some  local 
problems  of  size  (H/h)"^.  If /i  is  small,  we  cannot  have  both  l/H  and  H/h  small.  Thus  at 
least  one  of  the  subproblems  is  large.  The  computation  can  be  made  cheaper  by  recursively 
using  the  additive  Schwarz  method  to  solve  the  coarser  problems. 

We  consider  a  class  of  multilevel  algorithms  including  the  multilevel  additive  Schwarz 
method  considered  in  Dryja  and  Widlund  [28]  and  the  BPX  algorithm  studied  in  Bram- 
ble, Pasciak  and  Xu  [13]  and  Xu  [50].  We  use  the  space  decomposition  and  projection 
techniques,  introduced  by  Dryja  and  Widlund,  to  study  the  algorithms.  We  improve  the 
old  results  by  showing  that  the  condition  number  of  the  iteration  operator  grows  at  most 
linearly  with  the  number  of  levels  in  general,  and  is  bounded  by  a  constant  independent 
of  mesh  sizes  and  the  number  of  levels  under  the  ^^-regularity  assumption  of  the  elliptic 
operator.  We  note  that  similar  results  are  already  known  for  multigrid  methods. 

3.3.1      Description  of  the  Multilevel  Additive  Schwarz  Methods 

We  discuss  a  class  of  L-level  additive  Schwarz  methods.  We  define  a  sequence  of  nested 
triangulations  {7J=j}.  We  start  with  a  coarse  triangulation  T'  =  {t-,^},-'!,  where  t^  repre- 
sent an  individual  triangle.  The  successively  finer  triangulations  T'  =  {t'}{1  =  2,---,L) 
are  defined  by  dividing  each  triangle  in  the  triangulation  T'~'  into  several  triangles,  i.e. 

7-1  _  f   i)/Vi    rcfineqient    j2  _  l^^^N:,    refinement  refinement    jL  _  (^L-,Ni 
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We  assume  that  all  the  trlangulations  are  shape  regular.  Let  h\  =  diam(r,'),  hi  =  max,{/i'} 
and  h  =  h^. 

Let  V',l  =  !,•••,/,,  be  the  space  of  continuous  piecewise  linear  element  associated 
with  the  triangulation  T'.  The  finite  element  solution,  u/,  =  PyhU  £  V^  =  V^,  satisfies 

a{uH,<i>h)  =  f{4>h),     ^<f>heV''  =  V^.  (3.10) 

We  assume  that  there  are  L  —  1  sets  of  overlapping  subdomains  {17'},Ji,/  =  2,3,  ■•,L. 
On  each  level,  we  have  an  overlapping  decomposition 

We  assume  that  the  sets  {Cl\}  satisfy 

Assumption  3.3.1    The  decomposition  il  =  Uj_!ifi'  satisfies 

•  dCl\  aligns  with  the  boundaries  of  level  I  triangles,  i.e.  Cl\  is  the  union  of  level  I 
triangles.  Diameter{(l\)  =  0(/i/_i). 

•  On  each  level,  the  subdomains  {n'},Ji  form  a  finite  covering  of  ft,  with  a  covering 
constant  Nc,  i.e.  we  can  color  {J^'},-!i,  using  at  most  Nc  colors  in  such  a  way  that 
subdomains  of  the  same  color  are  disjoint. 

•  On  each  level,  associated  with  {n'},-Ji ,  there  exists  a  partition  of  unity  {tf'}  satisfying 

^$1  =  1,  with  el  e  Hoid'i)  n  C*(fii),0  <  e',  <  l  and  \V0',\  <  C/h,.i. 

i 

The  first  property  is  very  natural;  it  simply  says  that  the  restriction  of  the  triangulation 
T'  to  subdomain  Cl[  defines  a  triangulation  for  fi'  and  the  finite  element  problem  on  fi'  is 
well  defined.  The  second  condition  is  used  to  establish  the  upper  bound  of  the  spectrum 
of  the  additive  Schwarz  operator.  The  last  condition  is  used  for  the  lower  bound  of  the 
spectrum. 

One  way  of  constructing  subdomains  {S^'}, Jj ,  /  =  2,  •  •  • ,  L,  with  the  above  properties 
is  described  in  Dryja  and  Widlund  [24,25].  Each  triangle  r,'"* ,  t  =  1,  •  •  • ,  A^/,  /  =  2,  •  •  • ,  I, 
is  extended  to  a  larger  region  f/"^  so  that  ch'~^  <  dist(5f/"\5r/~^)  <  Ch'~^,  aligning 
5f,'~'  with  the  boundaries  of  level  /  triangles.  We  cut  ofi"  the  part  of  f'-~^  that  is  outside 
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fi.    We  use  f,'~'  as  the  subdomains  Q'-.    Another  way  of  constructing  {Q[}  is  given  in 
section  3.3.3. 

Let  iVi  =  1,  r,^  =  V^  and  V,'  =  V'  D  ^J(nJ)  for  i  =  1,  •  •  • ,  A',,  /  =  2.  ■  •• ,  i.  The  finite 
element  space  V''  =  V^  is  represented  by 

v'=Ev=i:Y:y,'-  (311) 

(=1    1=1  1=1 

Let  P'  :  V''  — ►  V-',  be  the  projection  defined  by 

aiPlu,<t>)  =  a{u,4>),     ^4>eV!. 

The  Z,-level  additive  Schwarz  operator  P  is  defined  by 

L    N, 
P  =  EEpI-  (3-12) 

/=0  1=1 

Instead  of  solving  the  original  finite  element  equation  (3.10),  we  solve  the  following  equiv- 
alent equation: 

Algorithm  3.3.1   Find  Uh.  €  V^  by  solving  iteraiively  the  equation 

Puh  =  9h 

with  gh  =  Yli  Jli9%-  f^^^  ^^^  9i  *"*  ^^^  solutions  for  the  following  finite  element  problems 

a{gl4>)  =  a{PyiuA)  =  f{4>l     ^AeV,'.  (3.13) 

I 

To  find  Uh,  we  first  find  the  right  hand  side  gh,  by  solving  (3.13),  and  we  then  use  the 
conjugate  gradient  method  to  sol^e  the  system.  In  each  iteration,  we  need  to  compute 
P/r/i  for  a  given  Vh,  G  V*  by  solving  the  equation 

aiPviVh,4>h)  =  aivh,4>h)    V<?1>/,  G  V/. 

This  is  a  finite  element  equation  on  Q[  with  mesh  size  hi,  and  dim(V/)  ss  c(hi_i/hi)^. 
Thus  the  size  of  all  such  problems  can  be  very  small. 


30 


3.3.2      Condition  Number  Estimate 

When  using  the  conjugate  gradient  method  to  solve  a  linear  system,  the  crucial  issue  is 
the  condition  number  of  the  iteration  operator.  Dryja  and  Widlund  [28]  established  the 
following  estimates  for  the  spectrum  of  P,  cf.  theorem  3.2  in  Dryja  and  Widlund  [28]. 

Theorem  3.3.1   For  P  defined  above,  the  following  inequalities  hold 

CiL-'^a(uh,Uh)  <  a{Puh,Uh)  <  C2La{uH,UH)  Vua  6  V''.  (3.14) 

Thus  k{P)  <  C^C^^L^,  i.e.  the  condition  number  of  P  grows  at  most  quadratically  with 
the  number  of  levels.  Here  the  constants  Ci  and  C2  are  independent  of  the  mesh  sizes  {hi} 
and  L. 

In  this  section,  we  improve  the  upper  bound  in  (3.14)  by  eliminating  the  dependence 
on  L.  Using  ^^-regularity,  which  holds  for  convex  regions,  we  can  also  eliminate  the 
dependence  on  L  in  the  lower  bound. 

Theorem  3.3.2  For  the  multilevel  additive  Schwarz  operator  P  defined  above,  we  have 

CxL~^a{uh,Uh)  <  a(Puk,Uh)  <  C-ia{uh,Uh)  Vu/i  6  V*" . 
If  the  equation  has  H^ -regularity,  then  the  lower  bound  can  also  be  improved,  and  we  have 

C\a{uh,Uh)  <  a{Puh,Uh)  <  C2a{uh,Uh)  Vu/,  6  V^. 
All  the  constants  are  independent  of  {hi}  and  L. 

Let  P'  :  V*"  -*'V',  be  the  orthogonal  projection.  Then,  for  u/,  £  V^  and  1  <  /  <  i, 

we  have 

/ 

P'uh  =  ^u',  with  u'  =  (P*  -  P'-^)uH,    u^  =  P^uh.  (3.15) 

1=1 

In  particular,  u/,  =  P^Uh  =  J2h=i  ^'-   Using  the  fact  that  V''  C  V',  for  A:  <  /,  we  obtain 
pkpi  ^  pipk  ^  pk    -pjjug  pi  _  pi-i  ig  ajso  a  projection,  and 

(P'  -  P'-')(P*=  -  P*-^)  =  %(P'  -  P'-^). 
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Therefore  ua  =  Yli'=o  "'  '^  an  a(-,  •)-orthogonaJ  decomposition  of  Uh,  and 

L 
a(«A,«A)  =  ^a(u',u')  (3.16) 

/=i 

Since,  on  each  level,  the  subdomains  {fi'LJi  form  a  finite  covering  of  fi  with  a  covering 
constant  Nc,  we  can  color  the  Q[  using  only  Nc  colors,  in  such  a  way  that  all  subdomains 
of  the  same  color  are  disjoint.  On  each  level,  we  can  group  the  subdomains  fi'  by  color, 
obtaining  Nc  sets  of  subregions.  Let  A,,  s  =  1,  •  •  ■ ,  Ad  be  the  sets  of  indices  of  these  sets. 
Since  the  subdomains  Cl[  of  the  same  color  are  disjoint,  subspaces  V^'  of  the  same  color  are 
mutually  orthogonal  and  Pj  P'  =  0  for  ii,t2  in  the  same  A,.  This  in  turn  implies  that 
^A.  =  E.6A.  P'  is  a  projection  from  V^''  onto  V^  =  E.ga,  V'-  ^  particular,  (P'^f  =  P^,. 
We  can  therefore  write 

1  =  1  3=]    l6A,  3-1 

where  the  Pj^    are  projections. 

The  next  lemma,  similar  to  the  strengthened  Cauchy  inequality,  plays  an  important 
role  in  obtaining  an  upper  bound  for  Xjn^{P). 

Lemma  3.3.1   For  k  <  I,  and  u*  G  V*,  we  have 

with  Tk,k  <  l,rk,k+i  ^  ^  ^'^^  ^k,i  <  Cr'-'^-'',  for  k  <  I  -  \. 

Proof.  Since  Pj^  is  a  projection,  the  conclusions  for  /  =  fc  and  I  =  k+l  are  trivial.  For 
/  -  1  >  /:,  we  decompose  A,  into  two  disjoint  sets  A;  and  Ag;  t  6  A/  if  the  subdomain  (l[ 
lies  in  the  interior  of  a  triangle  r*,  while  i  G  Ag  if  Cl[  intersects  drfi  for  some  t'.  We  write 
Pf^^  =  P\i  +  Psg.  We  note  that  u*  is  linear  in  each  triangle  r*  and  therefore  harmonic  in 
rj=.  For  each  «  G  A/,  P;'u*  6  ^o(f^')  C  HI{t^,)  for  some  i',  and  therefore,  a{P[y,  u*)  =  0. 
Thus, 

a{P'^u'y)  =  a{P[^u\u'') 

Let  S  =  supp{FAflU*}  =  U.eAsfi'-  Then 

aiP^y^u')  =  as(P'f,y,u'')  <  as(«*,u*) 
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Since  u*  is  linear  in  r*,  Vu*^  is  constant  in  r*.  Therefore 

,      ,        mes(5nr*)  ,,  /i,  ,  ,,  „  ,, 

J  mes(r*)        j  hk      '  j 

Summing  over  _;',  we  obtain 

Therefore 

I 
Proof  of  Theorem  3.3.2.  We  first  establish  the  upper  bound.  Since  V^^  =  J2ie\,  ^i  ^ 
V',  we  have  PJ^^P'  =  Pj,^.  Thus, 

aiP^,u,,Uf,)  =  a(P|.«;„Pj^,«,)  =  a{PiP'uH,PiP'uH). 

Substituting  (3.15)  into  the  above  equation,  we  obtain 

Jt=i  j=i  fc,i=i 

=  ii:\p'Ay\)'<c{j:n,\u''uf 

k=l  k=l 

Summing  over  all  colors  {I  <  s  <  Nc),  we  get 

Summing  over  /  and  changing  the  order  of  the  summation  for  k  and  /,  we  get 
L    Nt  t)        L     L 


^      '^   fc=i  l=k 


2  ^ 


<     CA'e(r^)^E  «(«''«*) 

=    C7Ve(i fa{uh,uh). 

1  —  r 
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In  the  last  step,  we  have  used  the  orthogonality  property  (3.16)  of  the  u^.  This  concludes 
the  proof  of  the  upper  bound  of  P  =  J2,i  Hi  P'- 

We  now  establish  the  lower  bound.  We  note  that,  for  the  general  case,  it  is  given  in 
Dryja  and  Widlund  [28].  When  the  problem  is  ^^-regular,  we  can  use  Nitsche's  trick  to 
show  that  the  a(-,  •)-projection  P'  satisfies  the  approximation  property 

.        ||/^«-t:|L.(n,  <C/i/|Hw(n)     Vu6^'(n).  (3.17) 

We  use  the  orthogonal  decomposition 

L 

U,  =  P'^UK  =  Y.u'  =  P\,  +  {P'-P')u,  +  .--  +  (P^  -  P^-'  )UH. 
1=1 

Since  u'  =  {P'  -  P'-')u/,  =  {P'  -  P'-^)^Uh  =  (/  -  P'-')u',  we  get,  using  (3.17), 

y\\mn)  <  C h,.,\u'\H^ti)-  (3.18) 

We  further  decompose  u'  as 

AT, 

u'  =  ^u; ,  with  «|  =  n'"(^|u')  e  v;'. 

Here  II'''  is  the  interpolation  operator  from  V^  to  V'^'  and  {$[}  a  partition  of  unity  as  in 
Assumption  3.3.1.  It  can  be  shown  that 

<    Ci\u'\l,^^,^  +  {l/hl,)\\u[\\l,^^,^). 

Summing  the  above  inequality  over  i,  using  the  finite  covering  property  of  {Cl[}  and 
inequality  (3.18),  we  obtain 

I  i  I 

<     C{\u%.^^^  +  l/hl,\\n'\\hm)  <  C\u%.^ay 

Summing  over  I,  I  <  I  <  L,  and  using  the  orthogonality  of  u',  we  get 

L 

IZZ]l"!l//'{n)  ^C'l"''l//>(fi)- 
/=i    1 

The  lower  bound  for  P  now  follows  from  Lions'  lemma.  I 
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Remark  3.3.1  Although  the  proof  of  the  upper  bound  is  given  for  the  model  problem, 
it  is  easy  to  see  that  it  works  for  any  uniform  elliptic  operator.  Since  we  can  confine  our 
study  to  one  substructure  at  a  time,  we  also  see  that  the  upper  bound  is  independent  of 
jumps  in  the  coefficients  between  the  substructures. 

3.3.3     A  Multilevel  Diagonal  Scaling 

We  begin  this  section  by  constructing  a  special  decomposition  of  the  domain  fi.  We  then 
show  that  this  decomposition,  and  the  corresponding  decomposition  of  the  finite  element 
subspaces,  satisfies  Assumption  3.3.1.  We  then  demonstrate  that  the  algorithm  is  a  mul- 
tilevel diagonal  scaling  (MDS),  a  natural  generalization  of  diagonal  scaling.  For  problems 
with  constant  coefficients  and  with  uniform  triangulations,  the  multilevel  diagonal  scal- 
ing algorithm  is  identical,  up  to  a  constant  multiple,  to  the  BPX  algorithm  of  Bramble, 
Pasciak  and  Xu  [13].  In  the  general  case,  BPX  with  diagonal  scaling  results  in  MDS 
algorithm.  ' 

Let  {T'}^j  be  a  nested  sequence  of  triangulations,  with  T'"*"^  obtained  from  T'  by 
dividing  the  triangles  (rectangles)  of  T'  into  four  triangles  (rectangles).  In  three  dimension, 
we  make  a  similar  construction.  We  consider  the  piecewise  linear  and  bilinear  elements 
or  trilinear  elements,  respectively.  As  in  the  previous  section,  the  finite  element  space 
associated  with  T'  is  denoted  by  V',  and  V'^  =  V^.  Let  4>\  be  a  nodal  basis  function  of 
V',  and  associate  with  each  <p'-  a  subdomain  (l\  =  supp{<^'}.  We  choose  V/  =  span{^'}  = 
F'  n  ^o(^i)  ^^^  obtain  the  decomposition 

L    N, 

and  the  projections  Pyi  :  V^  — ►  V^'.  Using  P  =  Ylf^-i  H,  A''  ^^  define  an  additive 
Schwarz  algorithm 

Algorithm  3.3.2  (MDS)  Find  the  finite  element  solution  u/,  6  V^''  by  solving  iteratively 
the  equation 

Puh  =  9h 

with  an  appropriate  right  hand  side  g^. 
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We  define  the  degree  of  a  vertex  i,  as  the  number  of  edges  directly  connected  to  x,, 
and  the  degree  of  a  triangulation  T**  as  the  maximum  of  the  degrees  of  its  vertices.  It  is 
easy  to  see  that  the  overlapping  subdomains  {Cl',}  satisfy  Assumption  3.3.1.  In  particular, 
we  see  that  on  each  level,  {fi'},^,  form  an  finite  covering  of  fi  with  a  covering  constant  less 
than  or  equal  to  degree(T')+l.  We  also  see  that  on  each  level.  {Cl[}  provides  relatively 
generous  overlap.  The  nodal  basis  {4>[]  can  be  chosen  as  a  partition  of  unity. 

The  optimal  convergence  properties  of  the  algorithm  follow  from  theorem  3.3.2,  and 
as  a  consequence,  we  obtain  an  improved  estimate  for  the  BPX  algorithm;  see  further 
discussion  below. 

To  see  that  the  above  algorithm  is  in  fact  a  generalization  of  the  diagonal  scaling 
method,  we  first  consider  a  matrix  representation  of  two  simple  algorithms. 

Algorithm  1.  In  the  two  level  additive  Schwarz  algorithm,  the  matrix  form  of  the 
projections  Pv,  is  given  by 

9 

after  a  permutation.  Here  A',  is  the  stiffness  matrix  associated  with  the  subspace  Vj.  Let 
n^  :  V^  —>  V'^  be  the  standard  interpolation  operator  and  Jlj^'  be  its  adjoint.  In  matrix 
form,  the  two  level  additive  Schwarz  algorithm  can  then  be  wTitten  as  a  preconditioned 
system 

where 

i=l 

Algorithm  2.  With  the  special  choice  of  the  subregions  fi,  =  supp{<^i},  and  Vi  — 
span{(^,),  the  additive  Schwarz  algorithm  corresponds  to 

where  D  =  diag(A'/,).  This  gives  us  the  Jacobi  conjugate  gradient  method. 

In  the  multilevel  case,  let  11,'  :  V''  — ►  V^'',(/i  >  /2)  be  the  standard  interpolation 
(prolongation)  operator,  and  let  n{|  :  V'^  -»  V'"  =  (Iljj )',(/)  >  /j)  be  a  local  averaging 
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operator,  which  is  the  adjoint  operator  of  11,'.  Algorithm  3.3.2  can  then  be  written  as: 
Find  the  solution  of  K^x  =  6  by  solving  the  preconditioned  system 

where 

Bi'  =  nf  A-f^ni  +  Ji^D^'nl  +  •  •  •  +  nt,Dil,ni-'  +  di\ 

Here  Ki  is  the  stiffness  matrix  associated  with  F'  and  Di  =  diag(A'/).  We  note  that  A'f ' 
can  be  replaced  by  any  good  preconditioner  B\  of  A'l. 

If  we  replace  the  matrices  Di  by  identity  matrices,  we  obtain  the  BPX  algorithm. 
However,  since  the  diagonal  elements  contain  information  on  the  shapes  of  the  triangles 
and  the  coefficients  of  the  problems,  we  expect  that  the  multilevel  diagonal  preconditioner 
will  work  better  in  practice  for  non-model  problems,  since  it  more  closely  reflects  the 
properties  of  the  problem. 

The  method  described  in  this  section  is  similar  to  the  hierarchical  basis  method  [52]. 
The  work  in  each  iteration  is  about  |  of  that  of  the  hierarchical  basis  method.  However,  the 
condition  number  is  much  better  than  for  the  hierarchical  basis  method,  and  the  method 
also  works  well  in  higher  dimensions,  at  least  for  problems  with  smooth  coefficients.  A 
detailed  comparison  of  the  BPX  algorithm  and  the  hierarchical  basis  method  is  given  in 
Yserentant  [53]. 

3.3.4     Numerical  Experiments 

In  this  section,  we  report  on  some  the  numerical  experiments  with  the  multilevel  additive 
Schwarz  methods. 

In  these  experiments,  we  only  concern  with  the  convergence  properties  of  the  algo- 
rithms. For  implementations  of  the  algorithms  on  a  parallel  computer,  see  Bjorstad, 
Moe  and  Skogen  [5]  and  Bj0rstad  and  Skogen  [6].  They  implemented  multilevel  additive 
Schwarz  algorithms  on  a  massively  parallel,  SIMD  machine  (MasPar  MP-1).  Approximate 
solvers  for  the  subproblems  are  also  discussed. 

The  experiments  were  carried  out  for  Poisson's  equation  on  a  unit  square  with  homo- 
geneous Dirichlet  boundary  conditions 


{ 


-Au    =     /(z)    in       fi,  ,        . 

1/    =       0       in    d^.  ^^-^^^ 
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total  no. 
of  elements 

no.  of  subdomains 
of  a  lower  level 

N  X  N 

ovlp  ratio 

no.  of  levels 
L 

cond  no. 

no.  of  iter, 
for  €  =  10-« 

8^ 

2x2 

1/2 

3 

7.2 

11 

16' 

2x2 

1/2 

4 

9.3 

17 

32'^ 

2x2 

1/2 

5 

10.7 

20 

64'^ 

2x2 

1/2 

6 

11.7 

21 

9^ 

3x3 

1/3 

2 

4.6 

9 

27'^ 

3x3 

1/3 

3 

7.1 

16 

81'^ 

3x3 

1/3 

4 

8.4 

19 

2-432 

3x3 

1/3 

5 

9.5 

21 

27' 

3x3 

1/3 

2 

4.8 

8 

81'^ 

3x3 

1/3 

2 

4.7 

7 

16^ 

4x4 

1/4 

2 

5.1 

13 

64^ 

4x4 

1/4 

3 

7.3 

17 

256^ 

4x4 

1/4 

4 

8.4 

20 

64'^ 

4x4 

1/4 

2 

5.3 

8 

25'^ 

5x5 

1/5 

2 

5.7 

14 

125'^ 

5x5 

1/5 

3 

7.6 

17 

Table  3.1:  Multilevel  Additive  Schwarz  Scheme,  Using  Bilinear  Element 

We  divide  the  domain  ft  into  N  x  N  square  elements  T}^,i,j  =  1,- ••,A'^,  and  obtain 
a  trianguJation  T'  =  {^A}-  We  then  divide  each  r/j  into  N  x  N  squares  to  obtain 
the  triangulation  7^  =  {r,^},  etc.  The  length  of  an  edge  of  t'^  is  denoted  by  Hi  and 
Hi  =  {l/N)'.  For  /  =  2,  •  •  • ,  L,  we  extend  r/~^  to  a  larger  square  f/~' .  The  overlap  ratio 

d\st{d  fi-\dri;'} 

Hi-i 


overlap  ratio  = 


measures  the  width  of  f/  '  \  r,'  '  in  terms  of  Hi^i,  the  side  of  the  square  t-j  ^.  We  use  Q 
as  our  domain  for  /  =  1,  and  ft',  =  f,~'  as  our  subdomains  for  /  =  2,  •  •  • ,  A'. 


In  these  experiments,  we  take  TV  =  2, 3, 4  or  5,  and  f-~^  \  r/~'  is  one  element  wide{.^;), 
i.e.  the  overlap  ratio  is  l/N .  Therefore,  we  only  need  to  solve  very  small  linear  systems 
of  order  9,  16  25  and  36,  respectively.  We  use  the  conjugate  gradient  method  to  solve  the 
system  Pu^  =  gn  iteratively.  The  last  column  of  the  table  gives  the  number  of  iterations 
required  to  decrease  the  /j  norm  of  the  residual  by  a  factor  of  f  =  10"^. 
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Chapter  4 

Some  Approximation  Properties 
of  Finite  Elements 


In  this  chapter,  we  discuss  some  properties  of  the  conforming  finite  element  for  the  bi- 
harmonic  problems.  We  first  discuss  and  prove  some  approximation  properties  of  the 
finite  elements.  We  then  define  some  discrete  norms  for  the  finite  element  spaces  and 
establish  their  equivalence  to  the  corresponding  continuous  Sobolev  norms  or  semi-norms. 
In  section  4.4,  we  estimate  the  norm  of  the  interpolant  of  a  product  of  two  functions. 
These  results  are  used  to  prove  the  estimates  for  the  convergence  properties  of  the  domain 
decomposition  methods  developed  in  chapters  5  and  6. 

Let  T''  =  {r,}  and  T^  =  {fi}  be  the  fine  and  coarse  triangulations  of  the  domain  fl, 
respectively.  Let  Vq,  V^  and  Vg  be  the  spaces  of  the  bicubic,  Argyris  and  Bell  elements, 
respectively,  associated  with  the  fine  triangulation.  Vq  ■,  V^  and  Vg  are  similarly  defined. 
For  a  vertex  i,,  let  r(a;,)  be  the  support  of  the  basis  functions  associated  with  i,,  i.e.  r(i,) 
is  a  region  containing  all  elements  with  this  vertex  Xi.  For  a  element  r,  let  f  be  the  smallest 
region  containing  r  and  all  its  neighboring  elements,  i.e.  f  =  U^  f-i^^^Tj. 

We  also  need  to  consider  some  subspaces  of  Vg,  V^  and  Vg.  Let  $f ,  \a\  =  0, 1,  and 
^^^  be  the  basis  functions  of  Vq  .  Let 

V(f  =  span{$f,|Q|  =  0,1,16  A''}. 

Thus,  Vq  is  the  subspace  of  Vq  consisting  of  functions  whose  mixed  second  derivatives 
vanish  at  all  vertices  of  the  substructures.  V^  C  V^  and  V^  C  V^  are  defined  similarly. 
We  use  h^{D)  and  A^{D)  to  denote  the  sets  of  fine  and  coarse  grid  points  in  the  set  D, 
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respectively. 

4.1      Approximation  by  Quasi-Interpolation 

In  this  section,  we  discuss  and  prove  a  stable  approximation  property  of  the  finite  element 
spaces  Vq  ,  V^  and  Vg.  The  results  are  based  only  on  the  two  properties  of  the  nodal  basis 
functions  stated  in  section  1.3.3.  Thus  if  another  basis  satisfies  the  same  assumptions,  the 
same  results  hold. 

We  describe  the  results  for  the  Bell  element;  the  results  for  the  .\rgyris  and  bicubic 
elements  are  similar. 

We  first  construct  a  local  linear  operator.  We  note  that  for  Q  C  7v^,  we  have  the 
inclusion  V^  C  V^  C  H'^(Q)  C  C^(n).  Let  <i>f,  \a\  <  2  and  let  r  be  an  element.  We  define 
a  local  linear  operator  Ig  :  H^{t)  — ►  Vb{t),  by 

/^uhX^u(i.)^.(^)  +  E  E{^""}^<^?(^)'      ^er.  (4.1) 

I  t     |a|=l 

Thus, 

(/^«)(i.)     =     u(x.), 
(Z?'*/^u)(x.)    =     {/?°u}.,  |q|  =  1, 
{D-rg){xi)    =     0,  |a|  =  2. 

We  find  that  if  p  is  a  linear  function  in  r,  then  Igp  =  p.  The  operator  Ig  has  the  following 
local  approximation  properties. 

Lemma  4.1.1   Operator  Ig  satisfies  the  following  inequalities 

|/b«Iw2(t)      <      C\u\„^r), 
l^«-«lw(r)       <      C/l^-'|u|//3(^),5  =  0,l. 

Ptxx)/.  By  Lemma  1.2.4,  we  know  that  there  exists  a  linear  function  p{x)  such  that  u 
can  be  written  as  p  +  /?  with  R  satisfying 

{D''R}r  =  0,|q|  =  0,1. 
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The  extended  Poincare's  inequality  gives 

\R\„.^r)  <  Crhl-'\R\„^r)  =  C./i^1u|h.(.).  (4.2) 

Using  the  linearity  of  Ig,  we  have  IgU  =  Igp  +  IgR  =  p  +  IgR.  By  the  definition  of  Jg 
and  using  the  fact  that  {D'^R}t  =  0  for  |q|  =  1,  we  get 

rgR  =  ^R{x,)4>,{x). 

By  the  triangle  inequality 

|/B^I//'(0<Ll^(^')ll<^>l//'(r)-  (4.3) 

From  Lemma  1.2.3,  with  k  =  l,\a\  =  0,il  =  t,  d  =  2,  and  (4.2),  we  have 

\Rix,)\  <  Crhr\R\„2(r)  =  C./i.|k|;^2(,). 

The  basis  functions  are  uniform  to  order  2  in  the  weak  sense,  i.e. 

\4>f\H'{r)  <  C/i^+'"'-%  for  \a\  =  0, 1,  \/3\  =  0,1,2, 

especially  |^,|;/i(t)  <  Ch}~' .  Substituting  the  above  two  inequalities  into  (4.3)  gives 

\lBR\H-(r)  <  ChMHHr)h\-'  <  C/l2-'|u|/,2(,),  5  =  0,1,2.  (4.4) 

The  first  part  of  the  lemma  follows  from  (4.4)  with  5  =  2.  By  the  triangle  inequality  and 
by  using  (4.2)  and  (4.4),  we  have 

IFgU  -  U\„.   =   IFgR  -  R\h.   <   irgRlH.  +  \R\h'   <  Ch'-'\u\„,^,y  (4.5) 

This  proves  the  second  part  of  the  lemma.  I 

We  now  construct  a  global  linear  operator.    Let  lyh  '■  H'^{il)  — ►  Vg,  be  the  quasi- 
interpolation  operator  defined  by 

i€A''(n)  ieA''(n)|c>|=i 


Thus, 


IyhU{x,)  =  Uh{Xi)  =  u{xi), 

D'^I^.uix,)  =  D^Uhix,)  =  {/?Mr(x.)>    l«l  =  1' 

D°Iy,u{x,)  =  D^'ukixi)  =  0,  |q|  =  2  , 

B 
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where  t{x,)  =  Ur.g-  '^j  defined  as  above.  The  quasi-interpolation  operator  Iy>,  :  H'^(Q)  — > 
V^,  and  /(^h  :  n^(n)  -  V^,  are  defined  similarly. 

The  quasi-interpolation  operator  lyh  has  the  following  approximation  properties. 

Theorem  4.1.1   Given  a  function  u  G  H'^{Q),  u/,  =  lynu  6  V^  satisfies 

l«/.|//2(n)     <    C\u\fJ7^Q), 
l"A  -  «lH'(n)     <    C/i^~'|u|;/2(n),5  =  0,1.2. 

^ere  the  constant  C  is  independent  of  h. 

Proof.  We  first  consider  one  element  r.  We  note  that,  locallv  in  an  element  r,  /,>h  is 

B 

defined  similarly  as  Ig,  except  that  1^,^  uses  the  average  of  the  first  deri\-atives  over  7-(i,), 

B 

while  Ig  uses  the  average  of  the  first  derivatives  over  r.  Let 


ef  =  {Z?'^ti}.(..)  -  {i?"u}„      |a|  =  1. 


Then 


|a|=l.eA''(T) 

<  E  E  kri  i<^r(x)iHV).  (4.6) 

|o|=lt€A''(T) 
From  the  fact  that  the  basis  functions  4>°{x)  are  uniform  to  order  2  in  the  weak  sense,  we 
have 

\4>f{x)\H.ir)  <  Cy^^'-\  for  \a\  =  1,5  =  1,2. 

Using  Poincar^'s  inequality,  we  obtain 

<    CA'-''/2|D'*«|;,,(,(,.))  <  Ch'-'lMm(,). 
Substituting  this  into  (4.6),  we  get 

|/(,.«  -  /^«|;,.(,)   <  Ch'-'\U\„,(,)       S  =  1.2. 

Combining  this  inequality  and  (4.5)  and  using  the  triangle  inequality,  we  have 

\Iy^U  -  U\„.^r)   <  \I^hU  -  Ilu\,f.^r)  +  \u  -  /fiUlH-lr)   <  C/l'~*|u|//3(f )  (4.7) 
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and 

\Iv^u\„2^r)   <   Uv^^  -  ^B^IhHt)  +  l^B^lffHr)   <  C\u\ff^fy  (4.8) 

B  B 

The  theorem  follows  by  summing  (4.7)  and  (4.8)  over  all  the  triangles  t  C  0.  and  using 
the  finite  covering  property  of  {f}.  I 

Note  that  in  Theorem  4.1.1,  we  made  no  conclusion  about  the  boundary  values  of 
lyhU.  However,  if  u  G  ^o(^)>  ^^  have  for  a  boundary  triangle  r  (r  n  90  /  0), 

|{«"}r-0|<C/l|u|;,.(n),   |a|  =  l. 

Using  this  fact  and  a  simple  modification  at  the  boundary,  we  can  show  the  following 

Theorem  4.1.2   Given  a  function  u  G  Hq{Q),  we  can  find  u^  6  V'^  PI  Hq{Q)  such  that 

WhlH^iO)     <     C|u|//2(n), 
|w/.  -  «l//'(n)     <    C/i^"''|u|H2(n),5  =  0,1,2. 

The  conclusions  of  Theorem  4.1.1  and  Theorem  4.1.2  also  hold  for  the  Argyris  and  bicubic 
elements,  if  we  replace  the  subspace  Vg  by  V^  and  Vq,  and  the  smooth  interpolation 
operator  /v>h  by  /,-,    and  /,-, 

*^B  •'/»  ''Q 

4.2      Approximation  by  Interpolation 

A  proof  of  the  following  lemma  can  be  found  in  Yserentant  [52]. 
Lemma  4.2.1 

-^  /  \uix)\dx  <  C{\og{Rla)  +  1/4)^/2,1^11^^^^^^ 

where  S  denotes  the  disk  of  radius  R>  a  >  0  with  center  xq  and  u  6  Hq{S). 

m 

We  know  that  for  linear  element 

Wh\l^(D)  <  0(1+ login /h))'f^\\u,\\H.^o), 

and  that  the  right  hand  side  can  be  replaced  by  the  semi  norm  |«a|//i(£))  if  the  mean  value 
of  u  over  D  vanishes;  cf.  Bramble  [9],  Bramble,  Pasciak  and  Schatz  [10]  and  Yserentant 
[52].  Similar  results  for  the  bicubic  element  V^^,  the  Argyris  element  V^  and  the  Bell 
element  Vg  can  be  established  by  applying  Lemma  4.2.1  to  Vu. 
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Lemma  4.2.2  For  u  €  V^,  we  have 

Hw..oo(n)  <  C(l  +  log(7///i))'/^||u,||„,(n). 

The  right  hand  side  can  be  replaced  by  the  semi-norm  \uh\nnD)>  '/{"}£>   -   {"x}d  = 
{"v}d  =  0. 

Let  V''  be  one  of  the  spaces  V^,V^  or  V^,  and  let  V^  denote  the  relevant  subspace 
Vq  ,  V^  or  Vq  .  Let  UyH  and  X[yn  be  the  standard  interpolation  operators  from  V^''  to 
V"  and  V",  defined  by 

|o|<2.eA«(n) 


and 

dx 


n<.HUA  =  E   E  ^«'.(^.)^r(^)'  ^  e  fi. 


|a|<li6A«(n) 

Lemma  4.2.3   The  interpolation  operator  IIv'h  has  the  following  properties: 
Ifuh  e  V^,  and  uh  =  UyHUh  6  V",  </jen 

l«A  -  nH\]f.^u)  <  cn'^'-'HH/h)^\u,\'„,^n)^  s  =  0, 1,2, 

and 

\uh  -  ^H\w..o.(n)  <  c^'<>-'H/r//i)'|«A|«j(fi),  5  =  0,1 

Proof.  Wegiveaproof  for  the  Bell  element.  Let  fi;  be  a  substructure,  and  i,,  i  =  1,2,3, 
be  its  vertices.  Let  $"  be  the  basis  functions  for  the  Bell  element  at  the  vertex  x,.  We 
have 

.eA«(n,)|o|<2 
We  want  to  estimate  the  difference  «;,  -  TlyHUh-    Using  the  linear  function  invariance 
property  of  UyH,  we  can  assume,  without  loss  of  generality,  that  {Z)"u}n,  =  0,  |q|  =  0, 1. 
By  the  triangle  inequality,  we  have 

|«//lw..»(n,)    <       E      E  l^"«(x.)||$r(r)|w..cc(n,) 

i6A«(n,)|o|<2 

•eA«(n,)|o|<2 
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and 

|«//l//.(n,)    <        E      E  l^"ti(x.)||$r(^  H'in,) 

.€A«(n,)|a|<2 
t'6A«(n,)|a|<2 

Using  the  fact  that  the  basis  functions  are  uniform  to  order  2.  i.e. 

\^n=')\w..^m<Cn\-\-'   and    |*°(x)|H.(n,)<C^'"'^'~% 
we  get 

and 

Since  {D^uyn,  =  0,  |a|  =  0, 1,  we  find  that 

|u|loc    <    CF|u|;/j(n,), 

l«lw'.-(n,)    <    C(l  +  log(F//i))»/Vl/P,Q,)' 
I«lw2.~(n,)    <    C'/i~^l"l//2(n,)- 

The  first  inequality  follows  by  the  embedding  Theorem,  Poincaxe's  inequality,  and  a  scaling 
argument.  The  second  is  just  a  consequence  of  Lemma  4.2.2.  The  third  is  the  inverse 
inequality.  Using  these  inequalities  in  the  estimate  for  lu/zlwo^.n,)  and  |u//|//»(n,),  we  get 

<  CH'-'{nih)\u\H,(^^), 

and 

<  cn'-'{n/h)\u\„^^,y 

By  the  triangle  inequality,  and  the  vanishing  mean  values  property  of  u  we  have 
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and 

The  conclusion  for  the  Bell  element  follows  by  summing  over  all  fi/.  The  results  for  the 
bicubic  element  and  Argyris  element  can  be  proved  similarly.  I 

Remark  4.2.1  It  is  easy  to  construct  examples  to  demonstrate  that  the  bound  in  the 
above  theorem  is  sharp.  As  in  the  3-d  problems,  interpolations  for  the  Hermitian  elements 
also  have  a  very  poor  bound. 

To  get  better  bounds,  we  consider  the  interpolation  interpolation  operator  TlyH-  The 
approximation  properties  of  Uyn  is  summarized  in  the  next  lemma. 

Lemma  4.2.4   The  interpolation  operator  UyH,  defined  above,  has  the  following  proper- 
ties: 
For  UheV",  UH  =  n^wUfc  eV",  we  have 

I"/.  -  "Hl//.(n.)  ^  CH'^'-'\l  +  login /h))\u,\l,^^^y 

and 

Proof.  We  prove  the  result  for  the  Bell  element;  the  assertions  for  the  bicubic  and 
Argyris  elements  can  be  proved  similarly.  Let  fi/  be  a  substructure,  and  i,,  i  =  1,2,3,  be 
its  vertices.  Let  $f  be  the  basis  functions  for  the  Bell  element  at  the  vertex  Xi.  We  have 

.eA«(n,)l'»l<i. 
Because  of  the  linear  function  invariance  property  of  11^,  we  can  assume,  without  loss  of 
generality,  that  {D°u}q,  =  0,|q|  =  0,1.   Note  that  the  only  difference  between  UyH^h 
and  UyHUh  is  that  UyH^h  does  not  contain  the  terms  D''u{x,)i°,\a\  =  2.  Therefore  we 
can  repeat  the  proof  of  Lemma  4.2.3  to  obtain 

|nv'H«|w.°c(n,)    ^        zl      z2  \^\w\''i°'(n,)\^?\\i''^(iit) 

•€A"(n,)|a|<l 
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<  C^-'|«|ioc(n,)  +  C^'-'\u\^^■^.oo^a,) 

<  CH'-'{l  +  \og{^/h)f/'\u\^p^^^y 

and 

|nv'H«|//.(n,)     <        ^      Yl  l"lw'i<'i.<^(n,)l</'ri//'(n,) 
ieA«(n,)|a|<i 

Using  the  triangle  inequality  and  the  vanishing  mean  value  property  of  u,  we  have 
Wh -^vHUh\w'-°°{Q,)    ^    l"lw.°°(n,)  +  |nv'HUAlH'».«'(n,) 

and 

<    Cn'-''il^log{H/h))\u\l,^^^y 
The  inequalities  for  operators  Ylyn  and  IlyH  can  be  established  similarly.  I 

4.3     Some  Discrete  Norms 

In  this  section,  we  define  some  discrete  norms  for  the  finite  element  spaces  Vq,V^  and 
Vg,  and  show  that  these  discrete  norms  are  equivalent  to  the  corresponding  continuous 
Sobolev  norms.  These  discrete  norms  are  used  in  establishing  the  optimal  convergence 
property  of  the  additive  Schwarz  method  for  the  biharmonic  equation  when  these  finite 
elements  are  used.  We  note  that  similar  ideas  can  be  used  to  show  the  optimal  property 
of  the  additive  Schwarz  methods  for  other  conforming  Lagrangian  or  Hermitian  elements. 
We  considered  the  bicubic  element  Vq,  the  Arg>-ris  element  V^  and  the  Bell  element 
Vg.  We  define  the  discrete  norms  in  terms  of  the  degrees  of  freedom  of  the  finite  element 
spaces. 

Definition  4.3.1    The  following  are  discrete  norms  for  the  finite  element  spaces: 

(1)  For  u  G  Vg.    Let  t  be  a  rectangle  with  vertices  Xi,i  =  1,2,3,4.    The  discrete  norms 
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over  T  are  given  by 

dxdy 

^dxdy 


|al  =  l      1 


,#j  |a|=l     t 

1  |q|=1  i/i  ■ 


u>/»ere 

Zu(i.)  =  u{xi)  -  {u(xi)  +  (Vu)(ii)  •  (x,  -  xi)}. 

^2^  For  u  £  Vg.  Let  r  be  a  triangle  with  vertices  i,.    The  discrete  norms  over  t  are  given 

by 

•/>  |or|=l      •  |o|=2     t 

I  |o|=li5tj  |o|=2     • 

(^5^  For  u  e  V^.  Z/Ct  t  be  a  triangle  with  vertices  i,.   The  discrete  norms  over  t  are  given 

by 

Ko     =     W\l,o+h'\D^'iy,)\\ 
l«ll,     =     \u\l^,  +  h'\D-'{y,)\\ 

I 

u;Aere  y,  =  (ajj  +  ijt)/2  »s  the  midpoint  of  the  edge  x^  x^. 

The  square  of  the  global  norms  are  defined  by  summing  the  square  of  the  local  norms 
over  all  T  C  il- 

Let  P{t)  be  a  space  of  polynomials  defined  on  r.  We  say  that  two  functionals  /i  and 
/j  on  P{t)  are  equivalent,  /i  ~  /2,  if  there  exist  two  constants  C\  and  C2  independent  of 
u  in  F(t)  such  that 

C,/,(ti)</2(u)<C2/i(u).  (4.9) 
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If  furthermore,  C\  and  C2  are  independent  of  r,  then  we  say  /i  and  /2  are  uniformly 
equivalent.  We  have  the  following  theorem  on  uniform  norm  equivalency. 

Theorem  4.3.1   Let  u  be  the  vector  consisting  of  the  degrees  of  freedom  of  the  finite 
element  in  a  triangle  r.   Then, 
for  ueV^ 

\\u\\l^^^     =     iRf^u,u)^\u\lo 


for  uEVb 


for  ueV^ 


The  theorem  still  holds  if  we  replace  the  norms  on  the  left  hand  side  by  the  norms  over  ft 
and  the  right  hand  side  by  sums  over  all  vertices  x,  G  SI. 

Proof  For  a  fixed  r,  the  left  and  right  hand  sides  define  quadratic  forms,  which  have 
the  same  null  space.  Thus  the  two  forms  are  equivalent,  since  they  are  defined  on  a  finite 
dimensional  space.  However  the  equivalency  constants  Ci  and  C2  may  depend  on  the 
triangle  r.  We  claim  that  Ci  and  C2  can  be  chosen  so  that  they  depend  only  on  the 
minimum  angle  of  triangle  r. 

This  is  trivial  for  the  bicubic  elements,  since  it  can  be  embedded  in  an  affine  family. 
Using  a  reference  element  and  a  mapping,  it  is  therefore  easy  to  see  that  the  discrete 
norms  are  uniformly  equivalent  to  the  corresponding  continuous  norms.  Argyris  and  Bell 
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elements  are  not  affine  elements.  To  check  the  uniform  norm  equi\-alency  for  the  Bell 
element  and  Argyris  element,  we  note  that  similar  triangles  have  the  same  equivalency 
constants.  Therefore,  it  is  enough  to  establish  the  result  for  triangles  of  unit  size  with 
a  minimum  angle  bounded  from  below  by  ^^in-  Let  AABC  be  such  a  triangle,  with 
A  =  (0,0),  jB  =  (1,0).  Then  its  third  vertex  C  must  lie  in  a  compact  region  D{9min)-  We 
note  that  the  continuous  and  discrete  norms  define  two  quadratic  forms  of  the  degrees  of 
freedom: 

|u|^  =  (A'ci,!)     and     \u\l  =  {Kjx,!), 

and  that  Kc  and  Kd  depend  continuously  on  C.  The  eigen\-alues  A(A'c)  and  X{Kd)  are 
continuous  functions  of  the  coefficients  of  the  matrices,  and  therefore  continuous  functions 
of  C.  Let  Amin(A'c)  and  XminiKd)  be  the  minimum  nonzero  eigenvalues,  and  Am4x(A'c)  and 
Anuui(A'j)  be  the  maximum  eigenvalues.  Since  C  lies  in  the  compact  set  D{6i^a),  there 
exist  constants  Ci  and  C2,  depending  only  on  Z?(tfmin).  thus  only  on  ^min)  such  that 

Ci   <  Anun(A'c)  <  An,4x(A'c)  <  Cj, 

and 

C\  <  Aniin(A'd)  <  AmAx(A'd)  <  €2- 

Since  the  matrices  Kc  and  Kd  have  the  same  null  space,  we  have 

(C,/C2)(A>,x)  <  iKdX,x)  <  {C2/Ci)iKcX,x),        VC  e  0(6^). 

Thus  we  have  proved  the  norm  equivalency  for  the  class  of  triangles  of  unit  size  with 
minimum  angles  bounded  from  below  by  ^min-  The  powers  of /i  in  the  discrete  norms  can 
be  obtained  by  shrinking  AABC  to  a  similar  triangle  of  diameter  h. 

The  global  norm  equivalency  follows  by  summing  over  all  elements.  I 

4.4      Estimates   for  the  Interpolant   of  a  Product   of  Two 
Functions 

In  the  proofs  of  the  estimates  for  the  domjiin  decomposition  methods,  we  need  to  estimate 
the  norm  of  the  interpolant  of  the  product  of  two  functions,  i.e.    |n;h(^u)|//2,  where  6 
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is  a  smooth  partition  function  and  «  is  a  finite  element  function.   For  simplification,  we 
assume  that  ${x)  £  C°°{Q).  Let 

00  =  |«(a:)|i,co(,),      Oi  =  |d(x)|vvi,oo(,),      02  =  |e(x)|H'2,oc(,).  (4.10) 

Let  Tl^rh,  riy^h  and  Ilvh  be  the  standard  interpolation  operator  to  the  finite  element  spaces 

Q  A  B 

Vq,V^  and  V"^,  respectively.  We  then  have  the  following  estimates 

Lemma  4.4.1 

\llv^{Ou)\l,^,)  <  C  l^el\Uy.u\j,,^,^  +  Ql\Uy.u\jf,^^^  +  0^||n^.^.t;||i,(,)}  ,     Vu  G  v<5. 

(4.11) 

(4.12) 

(4.13) 
j4//  cons<an<s  are  independent  of  mesh  size  h. 

Proof.    We  prove  inequalities  (4.12)  and  (4.13).    The  proof  for  (4.11)  is  similar.    By 
Theorem  4.3.1,  using  the  equivalent  discrete  norm,  we  have 

=    Ti  +  Ta  +  ra,  (4.14) 

where 

Lu{xi)  =  u{xi)  -  (u(ii)  -  Vu(x,)  •  (x,  -  xi)). 

Using  the  identity 

L{6u){x,)  =  e{x,)Lu{xi)  +  u{x,)L9{xi)  -  (^(xi)  -  d(xi))(u(x.)  -  u(xi)), 
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we  have 

t 

<     ^  ^  \e{x,)Lu{x,)\'  +  \uix,)Leix,)\'  +  |(^(x.)  -  e(x,))iu{x,)  -  u(xi))|' 

I 

I  is  a  Taylor's  expansion  and  we  have  \LO{xi)\  <  ©2/1^-  Since 

we  get  an  estimate  for  Ti, 

T,  <  C{Ql\Uy,u\j^,^,^  +  Ql\\lly.u\\mr)  +  0?|n;..«|^,(,)}.  (4.15) 

Manipulating  the  terms  in  the  expression  of  T2,  we  have,  by  elementary  calculus 

T2      =      J2Y1    \9{x,)UaiXi)-i-u(x,)ea{Xi)-e{xj)Uc{Xj)-u{x,)e^ix,)\^ 
•#J  |o|=l 

=     L  E  {(""C^-)  -  ^c,{x,))9{x,)  +  (^(x.)  -  ^(x,))u,(x,) 

t#J  |o|=l 

+(u(x.)  -  u(x,))^<>(x.)  +  (^^(x.)  -  e^{xj))u{x,)y 


<  c{0g|nv,.tz|^,,  +  0?|nv,.«|^,(,)  +  0^||nv,.u||2^,(,j}.  (4.16) 

T3  can  be  estimated  as 

|o|  =  2  |oi+crj|  =  2 

{o\=2  |a|=l  I 

<  0ol^v^«|?/.(.)  +  0?l^v.^u|?/.(,)  +  0^||^^.^u||i,,,).  (4.17) 

Inequality  (4.12)  follows  from  the  inequalities  (4.14)-(4.17). 
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Toestablishinequality(4.13),weonly  need  to  estimate  53,  |Z)"'(0u)(y,)-£)"'(^u)(i,)|' 
and  h^\D^'(6u){yi)\^ .  We  estimate  the  first  one;  the  estimate  of  the  second  is  similar. 

i 

=  E  l^n.(y.)"(y.)  +  ^(y.)«n.(!/.)  -  ^n.(x.)"(l.)  +  e{Xi)Un,{x,)\^ 

=  E  l^n.(l/.)(«(y.)  -  «(^.))  +  (^n.(y.)  -  e„.(x.))w(x,) 
t 

+  (fl(j/.)  -  e(x,))«n.(j/.)  +  ^(a:.)(«n.(l/.)  -  «n.(x.))l' 

+  C  Yl  mV')  -  ^(X.))«n.(y.)l'  +  CE  l^(^.)(«n.(j'.)  -  «n.(^.))l' 

=  Ti  +  Ta  +  Ta  +  T4.  (4.18) 

Using  Theorem  4.3.1,  each  term  can  be  estimated  as  follows: 

i 

7^2     <     0^/i2^u(i.)2<Cei|nv.^.tz|i,(,j.  (4.19) 

i 

T3    <     e?/i2^|u„.(y.)l'<C0?|nv,. «!!,,(,).  (4.20) 

t 

T4     <     e2^|u„.(y.)-«n.(x.)P<C02|nvjj«|?,2(,).  (4.21) 

i 

Using  the  basis  representation  for  functions  in  Vs,  we  have 

j  J     |a|<2  j 

Recalling  that  5Zf  <^i(x)  =  1,  using  the  fact  that  u  £  V5  and  the  fact  that  the  basis 
functions  for  the  Argyris  element  are  uniform  to  order  2,  we  have 

l"(y.)  -  «(x.)p 

=  lE[«(x,)-«(x.)]<^,(y.)  +  E  E  «a(x,)<^,"(y.)  +  E"n,(j/,)<^;'(y.)l' 

j  j     |a|<2  J 

<cEK^;)-«(^.)I'  +  E  E  /i"°'Mx,)|'  +  E'''K(y.)l' 

J  J     |a|<2  j 
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Thus, 

T,<Ql^\u{y,)-u{x,)\'<CQl\u\l,^^y  (4.22) 

I 

Inequality  (4.13)  now  follows  from  (4.18)-(4.22).  I 
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Chapter  5 

Additive  Schwarz  Methods  for 
the  Biharmonic  Problem 


In  this  chapter,  we  study  additive  Schwarz  methods  for  the  biharmonic  problem  using  the 
abstract  framework  of  chapter  2.  We  use  the  bicubic  element  Vq,  the  Argyris  element  V^ 
and  the  Bell  element  Vg. 

Suppose  that  the  finite  element  space  V  can  be  written  as  a  sum  of  subspaces 

where  Vq  is  a  coarse  subspace,  and  Vi  subspaces  associated  with  subregions  fi,.  Instead  of 
solving  the  original  finite  element  equation,  the  additive  Schwarz  method  is  introduced  in 
terms  of  an  auxiliary  problem:  Find  u  6  V  by  solving  iteratively  the  equation 

P«A  =  (/Vo  +  Pv,  +  •  •  •  +  PvJUh  =  9h 

for  some  g^. 

The  naturtd  question  is  how  to  find  decompositions  of  V  and  "what  properties  of  the 
decomposition  give  optimal  algorithms. 

As  we  pointed  out  earlier,  the  coarse  problem  is  crucial  in  our  algorithm.  In  the 
second  order  case,  an  obvious  candidate  for  the  coarse  subspace  is  the  space  associated 
with  T^ .  However  for  the  biharmonic  case,  when  Argyris  and  Bell's  elements  are  used, 
the  coarse  finite  element  spaces  V^  and  Vg  are  not  subspaces  of  the  fine  finite  element 
spaces  V^  and  Vg.  Therefore,  new  coarse  subspaces  have  to  be  found.  The  discovery  of 
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the  appropriate  coarse  subspaces  results  directly  from  an  obserN-ation  of  a  weak  coupling 
property  between  the  different  degrees  of  freedom  of  the  elements  and  some  results  on 
discrete  Sobolev  norms.  This  observation  also  leads  to  some  simplified  algorithms,  which 
turns  out  to  be  useful  even  for  the  algorithms  using  the  bicubic  element. 

In  the  iterative  substructuring  case,  the  situation  is  even  worse.  Even  when  we  use  the 
bicubic  element  and  Vq  C  Vq  can  be  used  as  coarse  subspace,  the  direct  generalization  of 
some  algorithms  for  second  order  problems  results  in  algorithms  with  condition  numbers 
which  grow  at  least  like  {H /h)^.  Better  algorithms  are  obtained  by  adding  certain  vertex 
spaces  to  the  space  decomposition;  cf.  chapter  6. 

The  difficulty  of  proving  the  optimality  of  the  algorithms  depends  on  the  presence  of 
high  order  derivatives  in  the  definition  of  the  elements.  The  tools  that  work  for  second 
order  equation  and  linear  element  cannot  be  used  here. 

5.1      The  Bicubic  Element 

5.1.1      Basic  Algorithnis 

We  first  triangulate  the  domain  fi  into  nonoverlapping  rectangles  fl,,  i  =  1,---,A^,  to 
obtain  the  coarse  triangulation  T^  =  {n,}f.  Then  each  rectangle  fi,  is  further  divided 
into  smaller  rectangles  r,  to  obtain  the  fine  triangulation  T^  =  {t,}.  We  aJso  decompose 
il  into  overlapping  subdomains  fi,,  i  =  1,  •••,A'^.  We  assume  that  the  decomposition 
n  =  U^ifi,  satisfies  Assumption  3.3.1.  Let  A^{D}  and  A"{D)  be  the  sets  of  the  fine  and 
coarse  grid  points  in  the  set  D,  respectively. 

We  use  two  finite  elemient  spaces  Vq  =  Vq  and  Vq,  which  are  the  bicubic  element 
associated  with  the  triangulations  T^  and  T^,  respectively.  In  addition,  we  use  the 
subspaces  V;  =  V^{(l,)  =  V^  n  H^{Cl,). 

» 

Let  Pv,  :Vq  —^  Vi,  be  the  a(-,  •)-orthogonal  projection  and  let 

N 
•=0 

We  have  the  following  additive  Schwarz  algorithm 
Algorithm  5.1.1   Find  u/,  €  V''  such  that 

Pnh  =  9H,  (5.1) 
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with  gh  —  5Z,  9i  ond  the  gi  given  by  the  solutions  of 

ai9i,4>h)  =  aiPv,u,<i>H)  =  fi4>k),     V<^/.  €  K-  (5.2) 

To  find  u/i,  we  first  find  the  right  hand  side  g^  by  solving  equations  (5.2).  We  then  use 
the  conjugate  gradient  method  to  solve  the  system.  In  each  iteration,  we  need  to  compute 
Pvh  for  some  element  r/,  6  V^.  This  is  done  in  the  following  steps,  which  can  be  carried 
out  in  parallel. 

•  Compute  PoVh  by  solving  K^x  =  6.  This  is  a  finite  element  problem  in  the  subspace 
V^.  Here  K"  is  the  stiffness  matrix  and  dim(V;^)  =  4A',  where  N  =  |A^(n)|,  is 
the  number  of  interior  coarse  grid  points. 

•  Compute  PiVh  by  solving  K-^x  =  b.  This  is  a  finite  element  problem  in  the  subdomain 
Cli.  Here  Ki  is  the  stiffness  matrix  and  dim(Vi)  =  4n,,  where  n,  =  |A''(n,)|.  is  the 
number  of  fine  grid  points  inside  fi,. 

5.1.2      Simplified  Algorithms 

It  is  of  course  desirable  to  reduce  the  work  in  each  iteration  without  decreasing  the  rate 
of  convergence.  Let  (j)f  and  $f  be  the  nodal  basis  functions  of  V^  and  Vq  ,  respectively. 
To  describe  our  alternative  algorithm,  we  need  to  introduce  some  new  subspaces. 

•  The  vertex  space  Vx*  =  sp^n{(l>l^ (x)} ,  for  each  vertex  n-  We  note  that  for  <i>  £  Vx^, 
T{xk)  =  supp(</>)  is  a  polygon  of  diameter  0(h). 

•  Vi  =  span{(^^(i),  \a\  <  1,  A:  €  A''(n,)}.  Vi  is  a  subspace  of  Vj  consisting  of  functions 
with  vanishing  mixed  second  derivatives  (^5*-  )  at  all  the  vertices. 

•  Vq  =  span{$f,  \a\  <  l,t  6  A^(il)}.  V^  is  a  subspace  of  Vq^  consisting  of  functions 
with  vanishing  mixed  second  derivative  {q^q   )  at  all  the  vertices  of  the  substructures. 

Algorithm  5.1.2  Find  u/,  G  Vq  by  solving 

PC'Kh  =  {Py»  +  E  ^v.  +  E  ^^* )"'«  =  9h  (5.3) 

i 

with  the  appropriate  right  hand  side  g^. 

57 


jfceA'' 


The  computation  of  j/,  is  similar  to  Algorithm  5.1.1.  In  each  conjugate  gradient  step 
of  Algorithm  5.1.2,  we  need  to  solve  three  sets  of  problems. 

1.  Compute  PynVh  by  solving  k"x  =  b.    This  is  a  finite  element  problem  in  the 

Q 

subspace  K^^.  Note  that  K"  is  a  principal  minor  of  A'^,  and  dimfV'^^)  =  \6\m{V^). 

2.  Compute  Py  v^  by  solving  K^x  =  b.  This  is  a  finite  element  problem  in  the  subdo- 
main  fi,.  K^  is  a  principal  minor  of  A'f  and  dim(V,)  =  ^dim(V;). 

3.  Compute  /V,  v/»  by  solving  K^^x  =  6,  with  x  and  b  scalars,  where  A'^^  =  a{4il^ ,  (j^^") 
is  the  1  by  1  stiffness  matrix  corresponding  to  the  degree  of  freedom  associated  the 
mixed  second  derivative  (^^7  )  at  the  vertex  Xk-  These  problems  are  very  easy. 

5.2      The  Argyris  and  Bell  Elements 

We  next  consider  the  Bell  and  Argyris  elements.  First,  we  present  the  basic  algorithms. 
We  then  describe  the  computationally  more  efficient  algorithms. 

5.2.1      Basic  Algorithms 

The  triangulations  T^,T^  and  the  subregions  Cli  are  defined  as  in  the  the  second  order 
case.  We  first  divide  the  domain  Q.  into  non-overlapping  substructures  fl,,  t  =  1,2,  ••  -A^. 
All  the  substructures  fi,  are  further  subdivided  into  elements  rj.  Thus  we  have  two  levels 
of  triangulations,  namely  the  coarse  triangulation  T^  -  {Q,}  and  the  fine  triangulation 
T''  =  {tj}.  We  denote  by  /i,  and  Hi  the  diameters  of  element  r,  and  substructure  fi,, 
respectively,  and  let  h.  =  max,  /i,  and  E  =  max,  Hi.  We  assume  that  all  the  substructures 
and  elements  are  shape  regular  in  the  usual  sense.  We  also  extend  each  substructure  to 
a  larger  region  Cl,.    We  assume  that  the  distance  between  the  boundaries  dQi  and  d(l, 

9 

is  bounded  from  below  by  a  fixed  fraction  of  Hi  and  that  50,  does  not  cut  through  any 
element.  It  is  dear  that  {fi,}  satisfies  Assumption  3.3.1. 

Let  Vg  and  V^  be  the  space  of  Bell  elements  with  respect  to  T''  and  T^ ,  respectively. 
V?  and  Vj^  are  similarly  defined.  We  present  the  algorithms  for  the  Bell  element.  The 
algorithms  for  the  Argyris  element  are  similar. 

In  general,  the  second  derivatives  of  $  €  Vg   at  the  edge  nodes  i,  have  two  diflferent 

58 


values  except  at  the  vertices  of  the  substructures.  Therefore.  Vg  ^  Vg.  Thus,  a  new 
coarse  space  has  to  be  found.  An  easy  way  of  modifying  V^  to  achieve  this  goal  is  by 
replacing  the  basis  functions  of  Vg  ■  Note  that,  in  a  substructure  flj,  basis  functions  $  of 
Vq    can  be  represented  by  the  basis  functions  of  Vq\ 

i,€A''(n,)|o|<2  x,6A''(an_,)|Q|<2 

which  is  now  replaced  by 

*(^)=    E    E*-(^')<^^w+    E    E  ^-(^')<^?(^)- 

x.eA''(nj)|o|<2  t€A''(anj)|o|<2 

Note  that  the  second  derivatives  of  ^  vanish  at  the  nodes  on  ^fi,.  We  define  the  coarse 
space  Ug  as 

{/^  =  span{*Ma|  <  2,  i^k"). 

It  is  easy  to  see  that  we  have  the  inclusion  U^  C  Vg-  Let  Vo  =  U§,  V;  =  V^  n  H^(Cli). 
Then  we  obtain  a  space  decomposition 

N 

1=0 

and  the  corresponding  orthogonal  projections  Py^. 
Algorithm  5.2.1   Find  Uh  €  Vg  such  that 

PuK  =  (Pvo  +  /V,  +  •  •  •  +  Pvs )«/.  =  9h  (5.4) 

with  gh  =  ^0  J,  and  gi  given  by  the  solutions  of 

ai9i,  4>h)  =  a{Pv,u,  <ph)  -  f{<f>h),     '^<t>h  e  Vi 

Once  we  have  computed  the  right-hand  side  gn,  we  use  the  conjugate  gradient  method 
to  solve  equation  (5.4).  In  each  iteration,  we  need  to  compute  Pvt^  for  some  element 
Vfi  6  V^.  This  is  done  in  the  following  steps 

1.  Compute  PoVh  by  solving  K^ x  =  6.  This  is  a  finite  element  problem  in  the  subspace 
Ug.  Here  K^  =  Kuh  =■  {a('J',,  *_,)}  is  the  stiffness  matrix  associate  with  Ug,  using 
the  modified  basis  {*f }.  We  note  that  dim(Vo)  =  dim(f/^)  =  dim(V^)  as  6iV, 
where  N  =  |A^(ft)|  is  the  number  of  substructure  vertices. 
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2.  Compute  P.r/,  by  solving  K^x  =  6.  This  is  a  finite  element  problem  in  the  subdomain 
n,.  The  number  of  unknowns  «  6n,,  where  hi  =  |A''(Q,)|  is  the  number  of  nodes 
inside  fi,. 

Remark  5.2.1  Notice  that  to  solve  the  coarse  problem,  we  need  to  form  the  matrix 
K"  =  KyH  =  {o(*,,  *;)},  where  ♦,  and  *j  Jire  the  modified  basis  functions.  Although 
the  stiffness  matrix  K"  =  KyH  =  {a($,,4»j)}  can  be  computed  explicitly  in  terms  of  the 
coordinates  of  the  three  vertices  of  fi/,  the  computation  of  Kuh  for  the  modified  basis  is 
not  straightforward.  A  standard  way  is  to  use  numerical  integration.  An  alternative  is 
to  replace  KyH  by  KyH.  This  is  equivalent  to  using  an  inexact  solver  Pq  to  replace  Pq. 

IB  B 

It  can  be  shown  that  {a($,,  $_,)}  and  {a($,,$j)}  are  spectrally  equivalent.  Thus  we  still 
have  an  optimaJ  algorithm. 

5.2.2      Simplified  Algorithms 

The  simplified  algorithms  for  Bell  and  Arg>ris  element  are  quite  similar  to  those  for  the 
bicubic  element. 

To  describe  these  algorithms,  we  need  to  introduce  the  following  subspaces. 

•  The  vertex  space  V^.^  =  span{<^^(i),  |q|  =  2},  for  each  vertex  Xk- 


•  Vi  =  span{</>^(i),  |q|  <!,/:€  A  (Q,)}.   Note  that  Vi  is  a  subspace  of  Vj  consisting 
of  functions  with  vanishing  second  derivatives  at  all  the  vertices  of  the  elements. 

9  tj^  =  spanftfja]  <  l,i  G  h"{9,)].    U§  is  a  subspace  of  U^  consisting  of  the 
functions  with  vanishing  second  derivatives  at  all  the  vertices  of  the  substructures. 

•  0^  =  span{*f ,  |q|  <  1,  * -jC^:),  i  6  A"iil)}-  V"  is  a  subspace  of  U^  consisting  of 
functions  with  vanishing  second  derivatives  at  aJl  the  vertices  of  the  substructures. 

Algorithm  5.2.2  Find  u/,  €  V^  such  that 

P^'W  =  iPo»+T,Pi'.+  Yl  Pv^, )«''  =  9H  (5.5) 

with  the  appropriate  right  hand  side  g/,. 


60 


The  computation  of  gn  is  similar  to  Algorithm  5.2.1.  In  each  iteration  step  of  Algo- 
rithm 5.2.2,  we  need  to  solve  three  sets  of  problems: 

1.  Compute  PyHVh  by  solving  K^x  =  b.  This  is  a  finite  element  problem  in  the 
subspace  Ug  .  Note  that  K^  is  a  principal  minor  of  K^  and  d\m{Ug)  «  ^d\m(Ug). 

2.  Compute  Py  Vh  by  solving  K^x  =  b.  This  is  a  finite  element  problem  on  the  subdo- 
main  ft,.  Here  K-^  is  a  principal  minor  of  K^  and  dim(Vi)  w  ^dim(V;). 

3.  Compute  Py^  Vh  at  each  interior  vertex  xjt  by  solving  K^^x  =  b.  Here  x,b  £71^  and 
■ft'rt  =  {a(*?»  *f )}  is  a  3  by  3  matrix,  *^  and  '*f  are  the  modified  basis  functions 
associated  with  the  second  derivatives  at  the  vertex  ijt.  These  problems  are  very 
easy. 

5.3      Optimality  of  the  BcLsic  Algorithms 

Our  result,  concerning  the  optimality  of  the  basic  additive  Schwarz  algorithm,  is  contained 
in  the  following  theorem. 

Theorem  5.3.1  The  condition  numbers  of  the  iteration  operator  P  in  Algorithms  5.1.1 
and  5.2.1  are  bounded  by  a  constant  independent  of  number  of  substructures  and  the  mesh 
sizes.  More  precisely,  we  have 

Cia{u,u)  <  a{Pu,u)  <  C2a{u,u),     Vu  6  V'', 

Proof  The  upper  bound  follows  from  the  finite  covering  property  of  {17,}  and  the 
properties  of  the  projections,  with  Cj  =  {Nc  +  1).  Here  Nc  is  the  finite  covering  constant. 

To  establish  the  lower  bound,  we  need  to  construct  a  good  partition  u  =  $2,  «,,  Vu  £  V^ 
in  order  to  use  the  Lions'  Lemma. 

For  the  bicubic  elements,  by  Theorem  4.1.2,  we  know  that  'iu^  £  V^  C  H^iH),  we  can 
find  SLUM  £V^  (^  -^oC^)  satisfying 

\uh\hhq)  ^  C:\uh\H^n)^  (5-6) 

and 

l«H  -  «A|//'(n)  <  C^^~''\uh\H^n)^  (5-7) 
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We  define  Uo  =  «//,  w^  =  Uh  -  u//.   Since  V^  C  V^ ,  we  find  uo  £  V^  C  V^.  Because 
of  the  generous  overlap  of  (l„  we  know  that  there  exists  a  partition  of  unity  {6,}  with 

9,  e  C^(n,),  which  satisfies 

N 

and 

e.  =  max|<?.|j^.,«(ft,)  <  CH-',  s  <  2. 

Let  n''  =  UvH  be  the  standard  interpolation  operator  to  the  finite  element  space  VA.  We 
define  u^  by 

u.  =  U\9,wt,). 

It  is  easy  to  check  that  u,  is  well  defined,  and  that  u,  €  V,.   By  the  linearity  of  n'',  we 

know  that 

N 

t=0 

From  Lemma  4.4.1,  we  know 

Using  the  fact  that  u;;,  =  u^  -  «;//  €  V'',  we  obtain 

Summing  over  r  G  fi,,  we  get 

The  conclusion 

.=0 

follows  from  the  fact  that  0,  <  CH'',  |u^/,|w(n)  <  ^^"'^/.l/z^in)  3-"^  the  finite  covering 
property  of  J),. 

For  the  Argyris  and  Bell  elements,  we  note,  from  the  results  on  the  discrete  Sobolev 
norms,  and  the  construction  of  Ug  and  t/^,  that  the  modified  basis  functions  ^°  have 
similar  properties  as  the  basis  functions  of  V^ .  In  particular,  * f  have  the  two  properties 
of  studied  in  chapter  1.  Thus,  the  stable  approximation  properties  hold  for  U§ .  The  rest 
of  the  proof  is  similar  to  the  case  for  the  bicubic  element.  I 
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5.4     Optimality  of  the  Simplified  Algorithms 

We  have  the  following  result  on  the  optimality  of  the  simplified  algorithms. 

Theorem  5.4.1  Algorithms  5.1.2  and  5.2.2  are  optimal,  i.e.  the  iteration  operator  P<^) 
satisfies 

Cia{u,u)  <  a(p(^^u,u)  <  C2a{u,u),  Vu  £  V^, 

where  C\  and  C2  are  independent  of  the  number  of  substructures  and  the  number  of  ele- 
ments. 

Remark  5.4.1  The  optimality  of  Algorithms  5.1.2  and  5.2.2  can  essentially  be  proved 
by  showing  that  the  vertex  spaces  Vx^  are  weakly  coupled  to  the  other  subspaces.  This 
property  will  also  be  important  in  developing  iterative  substructuring  algorithms  for  the 
biharmonic  equation. 

Proof.  We  know  that  {J7,}  forms  a  finite  covering  of  Q  with  a  covering  constant  Nc. 
It  is  easy  to  see  that  {T{xk)}ke\h  forms  a  finite  covering  of  Q  with  a  covering  constant 
degree(T'')  +  1.  Thus  {Cli,T(xk)}i^k  «Jso  forms  a  finite  covering  of  Q  with  a  covering 
constant  (degree(T'')  +  1)  +  Nc.  The  upper  bound  of  P  follows  with  C2  =  (degre€(T'')  + 
1)  +  A'c  +  1- 

To  establish  the  lower  bound,  we  need  to  construct  a  good  decomposition  u;,  =  Yli  "t 
in  order  to  use  the  Lions'  Lemma.  By  Theorem  4.1.2,  we  can  find  uh  6  Vq  D  ^o(^) 
satisfying 

l«//l//2(n)  <  C\uh\fP(n), 
and 

I"//  -  w/.l/f'(n)  <  Cn^-'\uh\fp^Q). 

Let  uo  =  U//5  ^h  =  t^h  -  ^H  and 


«z.(x)  =  ^^^h^4>iy(x),  for  X,  e  A\Q). 
From  the  properties  of  the  basis  functions,  it  is  clear  that 

l«xjt(a:)|//'(n)  <  C\wh{x)\H-(r(x^))- 
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Summing  over  all  ijt  6  A**,  we  have 

Taking  s  =  2,  and  using  the  fact  that  \wh\H^  <  C|ua|//2,  we  obtain 

EKJ^<CKI?/.(n).  (5.8) 

Let  Wh  -  rvh  -  E*"x»-    Then  |%|//.  <  C\wh\H--    Let  w,(i)  =  UyH{6,Wh)-    From  the 
linearity  of  Uyh,  we  obtain 

and 

D^'iuH  -u„-  J2  «xJ(^t)  =  0.  for  l«l  =  2,fc  G  A\ 
fceA" 

Thus,  Vuiixk)  =  0,  for  k  G  A'',  which  implies  that  u,  G  V",.  Using  the  same  technique  as 
in  the  proof  of  Theorem  5.3.1,  we  can  prove  that 

YlMh<C\w,\l..  (5.9) 

t 

Combining  (5.6),  (5.8)  and  (5.9),  we  obtain 

The  theorem  now  follows  from  Lions'  Lemma.  The  conclusions  for  the  Argyris  and  Bell 
elements  can  be  established  similarly.  I 

5.5     Numerical  Results 

In  this  section,  we  report  on  some  numerical  results  with  the  additive  Schwarz  methods 
for  the  biharmonic  equation.  The  computations  were  carried  out  on  a  Sparc  Station  1  and 
a  CONVEX  C-1  machine.  All  the  experiments  are  for  the  unit  square  domain  Q  =  [0, 1]^. 
We  divide  the  domain  il  into  Nh  X  Nh  idential  square  subdomains  fl,j  and  obtain  a 
coarse  triangulation  T"  =  {n,j}.  The  length  of  sides  of  fi,j  is  H  =  -^  and 

Q„  =  {{i-\)H,iH)yi{j-\)HJH). 
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Each  square  ft,j  is  further  divided  into  smaller  square  elements  ti^  of  the  same  size  and 
we  obtain  the  fine  triangulation  T^  =  {r,j}.  Let  n\  be  the  number  of  elements.  Then  the 
length  of  sides  of  the  elements  is  h  =  -^  and 

Ti,  =  iii-l)h,ih)x{(j-l)h,jh). 

We  extend  each  fi.j  to  a  larger  region  ft,j. 

Clij     =     ({i-l)n-rH,iH  +  rn)xiiJ-l)n -rHJH  +  rH)  (5.10) 

=     {{i-l)n -meh,in  +  meh)x{{j-l)H -m.hJH  +  nieh).       (5.11) 

We  cut  off  any  part  of  Clij  that  is  outside  of  Cl.  Here  r  is  the  overlap  ratio  and  me  is  the 
number  of  elements  by  which  Qij  is  extended  in  each  direction.  We  note  that  rH  has  to 
be  a  multiple  of  h,  since  we  want  the  dClij  to  align  with  element  boundaries. 

For  simplicity,  we  want  all  the  subproblems  to  have  the  same  size,  i.e.  all  the  submaJns 
Clij  have  the  same  size.  For  this  to  hold,  we  need  to  modify  the  boundary  subdomains 
Clij.  Instead  of  extending  ilij  by  rrie  elements,  it  is  extended  by  2me  elements  towards  the 
interior.  This  is  especially  useful  for  computations  on  SIMD  machines;  cf.  Bjorstad  and 
Skogen  [6].  For  all  the  experiments  of  this  section,  we  stop  the  iteration  when  the  norm 
of  the  residual  is  reduced  by  a  factor  of  e  =  10"^. 

In  the  first  set  of  experiments,  we  discretize  the  equation  using  the  bicubic  elements 
Vq  and  we  use  Algorithm  5.1.1.  The  total  number  of  degrees  of  freedom,  dim(VQ)  = 
4{nh  -  ly.  The  results  are  summarized  in  table  5.1.  The  first  column  contains  4(n/,  —  1)^, 
the  total  number  of  degrees  of  freedom.  The  second  column  contains  Njf,  the  number  of 
subdomains.  In  the  third  column,  we  give  the  overlap  ratio  r  and  the  number  of  elements 
TTif.  by  which  (lij  is  extended  from  Qij.  We  note  that  the  real  overlap  ratio  for  boundary 
subdomains  will  be  slightly  different  from  that  of  the  interior  ones.  The  fourth,  fifth  and 
sixth  columns  give  the  minimum  and  minimum  eigenvalues  and  the  condition  number  of 
the  iteration  operator  P,  respectively.  The  last  column  contains  the  number  of  iteration 
required  to  reduce  the  norm  of  the  residual  by  a  factor  of  e  =  lO"**. 

We  observe  from  table  5.1  that,  in  general,  Xmax  is  between  four  and  five,  except  in 
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total  #  unkns 

#  of  subdom 

ovlp  ratio 

^min 

-^mAX 

K(P) 

#  of  iter. 

4  X  (8-  1)'^ 

22 

1,1/2 

0.99 

5.0 

5.0 

8 

4  X  (16-  1)'^ 

2^ 

1,1/4 

0.51 

4.3 

8.3 

10 

4  X  (16-  l)'' 

42 

1,1/4 

0.65 

4.5 

7.0 

11 

4x  (16-  1)'^ 

42 

2,1/2 

1.00 

9.0 

9.0 

12 

4x  (32-  1)^ 

42 

2,1/4 

0.65 

4.5 

6.9 

10 

4x(32-  1)'^ 

82 

1,1/4 

0.62 

4.5 

7.3 

11 

4x  (64-  1)'^ 

82 

2,1/4 

0.62 

4.5 

7.2 

10 

4x  (64-  1)^ 

162 

1,1/4 

0.61 

4.5 

7.3 

11 

4  X  (80-  1)2 

82 

2,1/5 

0.48 

4.3 

9.0 

11 

4  X  (80-  1)'^ 

16'^ 

1,1/5 

0.47 

4.3 

9.1 

12 

4  X  (100-  1)2 

102 

2,1/5 

0.48 

4.3 

9.0 

11 

4x  (100-  1)2 

202 

1,1/5 

0.47 

4.3 

9.1 

12 

Table  5.1:  ASM  Using  the  Bicubic  Element  V^. 

the  fourth  row,  where  A^ax  =  9-  The  reason  that  Xmux  =  9  in  the  fourth  row  is  that,  in 
this  case,  we  extend  the  n.j  too  much  and  the  covering  constant  Ncover  =  8. 

For  comparison,  we  also  compute  the  condition  numbers  of  the  additive  Schwarz 
method  without  a  coarse  subspace.  For  a  problem  with  16  by  16  grid  points,  4  by  4 
subregions  and  overlap  ratio  r  =  1/4,  we  have  k{P)  =  250. 

The  condition  numbers  of  the  stiffness  matrices  A'''  for  these  fourth  order  problems 
are  extremely  large.  We  give  two  examples  here.  For  n/j  =  10,  we  have  k(A'i2,i2)  *  1800 
and  for  n/,  =  20  we  have  k(A'2o,2o)  *  14000. 

In  the  next  set  of  experiments,  we  consider  Algorithm  5.1.2.  We  also  use  the  bicubic 
element.  The  two  triangulations  T'',  T"  and  the  subdomains  Cl,j  are  defined  as  in  the 
previous  case  and  so  are  the  overlap  ratios  r  and  m,.  The  operator  for  Algorithm  5.1.2  is 
defined  by 

The  Vq  is  the  bicubic  element  associated  with  the  coarse  triangulation  T",  Vj  the  bicubic 
elements  associated  a  subdomajn  ilki  and  V^^  the  vertex  spaces  associated  with  a  fine  grid 
point  Xki  =  (kh,lh). 
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total  #  unkns 

#  of  subdom 

ovlp  ratio 

'^mjn 

'^mAx 

k{P) 

4x(8-  1)'^ 

22 

1/2 

0.93 

4.49 

4.8 

4x(16-iy^ 

2^ 

1/2 

0.83 

4.47 

5.4 

4x(16-  ly 

2^ 

1/2 

0.65 

4.24 

6.4 

4x(32-  ly 

2^ 

1/4 

0.83 

4.45 

5.3 

4x(32-  ly^ 

4' 

1/4 

0.74 

4.43 

6.0 

4  X  (32-  1)'^ 

8^ 

1/4 

0.60 

4.43 

7.4 

4x(64-  ly 

8^ 

1/4 

0.64 

4.42 

6.9 

4x(64-  ly 

162 

1/4 

0.60 

4.41 

7.4 

Table  5.2:  ASM  Using  the  Bicubic  Element  V^:  The  Simplified  Version 

We  note  that  we  can  use  different  weights  for  different  projections.  In  our  experiments, 
we  use 

i 

The  results  are  summarized  in  table  5.2.  We  observe  that,  for  the  same  triangulations 
T''  and  T^ ,  the  condition  number  k{P^'^^)  is  approximately  of  the  same  size  as  k{P)  of 
Algorithm  5.1.1,  in  some  cases,  even  smaller. 

In  a  third  set  of  experiments,  we  consider  the  Bell  element  discretization.  We  use 
Algorithm  5.2.1  to  solve  the  linear  system.  Since  we  need  triangular  elements,  the  two 
triangulations  T^,  T^  are  defined  by  dividing  the  squares  obtained  from  the  previous 
cases  into  two  triangles.  For  simplicity  we  use  square  subdomains  n,j  as  in  the  previous 
cases.  The  results  are  summarized  in  table  5.3. 

We  recall  that,  in  this  case,  the  coarse  problem 

corresponds  to  the  linear  system 

K,thx  =  b. 

Here  V^  is  the  modified  coarse  space.  In  our  experiments,  we  do  not  form  matrix  ^v^- 
Instead,  we  use  K^h,  the  stiffness  matrix  for  the  coarse  Bell  elements.  This  is  equivalent 

B 

to  solving  Kijht  =  6  by  an  approximate  solver.    Since  K^h  and  K^h  are  spectrally 
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total  #  unkns 

#  of  subdom 

ovlp  ratio 

^mJn 

^m&x 

K(P) 

#  of  iter. 

6x  (16-  1)^ 

2^ 

2,1/2 

1.0 

6.0 

6.0 

9 

6x  (16-  1)2 

22 

1,1/4 

0.33 

4.3 

13.0 

12 

6x  (16-  ly 

4'^ 

1,1/4 

0.40 

4.6 

11.3 

14 

6x  (32-  1)2 

42 

2,1/4 

0.44 

4.5 

10.3 

12 

6  X  (32-  1)2 

82 

1,1/4 

0.35 

4.6 

12.8 

15 

6x  (64-  1)2 

8-2 

2,1/4 

0.38 

4.5 

11.9 

13 

6x  (64-  1)2 

16'^ 

1,1/4 

0.35 

4.6 

13.2 

14 

6x  (80-  1)2 

102 

2,1/5 

0.41 

5.4 

12.8 

13 

Table  5.3:  ASM  Using  the  Bell  Element  V^ 

equivalent,  we  still  have  an  algorithm  with  optimal  convergence  properties.  The  numerical 
results  are  summarized  in  table  5.3.  Since  we  do  not  solve  the  coarse  problem  exactly, 
^muxiP)  does  not  always  lies  in  between  4  and  5. 

5.6      Multilevel  Methods  for  the  Biharmonic  Problem 

In  this  section,  we  study  multilevel  additive  Schwarz  methods  for  the  biharmonic  equa- 
tion. Although,  all  the  multilevel  methods  studied  in  chapter  3  can  be  modified  for  the 
biharmonic  problem,  we  only  consider  a  special  case,  which  in  matrix  form  corresponds 
to  a  multilevel  block  diagonal  scaling  (MBDS).  To  simplify  the  presentation,  we  use  the 
bicubic  elements. 

We  first  define  a  sequence  of  nested  rectangular  triangulations  {T'};^j  as  in  case  of 
second  order  problems,  cf.  chapter  3.  We  use  r,'  to  denote  the  elements  in  T'.  The  level  / 
grid  points  are  denoted  by  A',  and  the  basis  function  by  4>[  ^,  i  6  A'. 

Let  V'  =  Vq'  be  the  bicubi«  element  associated  with  T'.  The  finite  element  solution 
UA  G  V^  satisfies 

a{uh,4>h)  =  KM,   <i>KeV''  =  V^.  (5.12) 

Let  il[  =  supp{<^'  „}  be  the  support  of  an  individual  basis  function  and  let  V-'  =  span{<^'  £,} 
be  the  span  of  the  level  /  basis  functions  at  the  grid  point  i^.  We  note  that  for  bicubic 
elements  dimfV/}   =  4.    On  each  level,  we  have  an  overlapping  decomposition  of  the 
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domain 

n  =  u,n|. 

This  decomposition  satisfies  Assumption  3.3.1.  In  particular,  the  decomposition  has  the 
finite  covering  property  and  there  exists  a  partition  of  unity  {8[}  associate  with  {il[}-  We 
use  the  space  decomposition 

L    N, 

1=1  t=i 
The  operator  of  the  L-level  additive  Schwarz  algorithm  is  given  by 

L    Si  L    N, 

i=\  «=i       /=i «=i 

with  Pyi  :  V'*  -*  F/,  the  a(-,  •)-orthogonal  projections. 
Algorithm  5.6.1  (MBDS  Algorithm)  Find  Uk  G  V^  by  solving 

PuH  =  9h  (5.13) 

with  an  appropriate  right  hand  side  gh  ■ 

In  matrix  form,  equation  (5.13)  can  be  written  as: 

where 

B-'  =  Di'  +  Hi,.,  ViHLi  +  •  •  •  +  niA'r^n'i. 

Here  Ki  is  the  stiffness  matrix  associated  with  T',  Di  =  diaglA"/),  11/  a  prolongation 
operator  and  II'  a  restriction  operator.  The  matrix  iff  ^  can  be  replaced  by  any  good 
preconditioner  of^A'i. 

Theorem  5.6.1    The  multilevel  additive  Schwarz  operator  P  satisfies 

CiL-^a(u, u)  <  a{Pu, u)  <  CLa{u, u). 

Thus 

k{B-'^K)  <  CL^ 
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Remark  5.6.1  For  second  order  problems,  the  corresponding  algorithm  is  a  multilevel 
diagonal  scaling  method,  which  is  equivalent  to  the  BPX  algorithm  of  Bramble,  Pasciak 
and  Xu  [13].  We  note  that  in  the  second  order  case,  we  obtain  a  constant  upper  bound  for 
the  operator  P\  cf.  chapter  3  and  Zhang  [54].  It  is  jJso  possible  to  strengthen  the  result 
of  Theorem  5.6.1. 
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Chapter  6 

Iterative  Substructuring 
Algorithms 

6.1      Introduction  and  Overview 

In  this  chapter,  we  study  iterative  substructuring  methods  for  the  biharmonic  problem. 
We  first  review  some  methods  for  second  order  problems,  using  the  additive  Schwarz 
framework  introduced  by  Dryja  and  Widlund  [25].  We  then  generalize  some  iterative 
substructuring  schemes  to  the  biharmonic  equation.  We  demonstrate  that  direct  general- 
izations of  the  iterative  substructuring  methods  designed  for  second  order  problems  give 
slow  algorithms.  Better  algorithms  are  obtained  by  adding  certain  vertex  spaces  to  the 
decomposition  of  the  finite  element  spaces. 

We  recall  that  iterative  substructuring  methods  are  domain  decomposition  methods 
using  non-overlapping  subregions;  cf.  Bramble,  Pasciak  and  Schatz  [10,11],  Widlund  [47]. 
Techniques  for  analyzing  some  iterative  substructuring  methods  in  the  additive  Schwarz 
framework  were  developed  by  Dryja  and  Widlund  [25,26],  where  it  was  shown  that  the 
algorithm  of  [10]  can  be  viewed  as  an  additive  Schwarz  method  with  a  set  of  specially 
chosen  subdomains. 

As  always,  we  use  two  triangulations  T^  =  {ft,}  and  T'*  =  {r,}.  Let  r,_,  =  x,j;_, 
be  a  common  edge  of  the  two  neighboring  substructures  fi,  and  ilj  and  let  ftij  =  ft,  U 
Tij  U  ftj.  The  subdomains  {ftij}  play  the  role  of  the  overlapping  subdomains  {ft,}  in  the 
decomposition  of  domain.  Using  {^ij},  we  obtain  a  decomposition  of  the  domain  ft. 

ft  =  Uftij. 
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Let  K'',  V^  be  the  linear  element  associated  with  T^  and  T^ ,  respectively.   Associated 
with  each  edge  r.j,  we  have  a  subspace  Vij,  defined  by 

Using  the  space  decomposition 

•J 
and  the  corresponding  projections 

we  obtain 

Algorithm  6.1.1   Find  Uh  €  V^  by  solving 

Puh  =  9h 

with  an  apprvpriate  right  hand  side  g^- 

Let  Vj-^*  C  Vij  be  space  of  discrete  harmonic  functions  and  let  V,  =  V^  d  /foCfii). 
Then  we  have  an  orthogonal  decomposition  of  Vij: 

K,  =  K  +  V,  +  v;)^. 

Thus,  iV.j  =  Pv,  +  Pvj  +  P  ha"   Computing  Pv,Vh  corresponds  to  solving  a  local  finite 

•J 
element  problem  in  fl,.  In  matrix  form,  P  jj^  is  closely  related  to  the  Schur  complement 

'} 
Sij.   Let  J  be  the  square  root  of  the  discrete  Laplacian  operator  along  the  edge  r,j.   If 

we  replace  Sij  by  J  in  our  preconditioner,  we  obtain  the  BPS  algorithm  considered  in 

Bramble,  Pasciak  and  Schatz  [10]. 

The  condition  number  of  the  BPS  algorithm  grows  like  (1  +  log^{H/h)),  see  Bramble, 
Pasciak  and  Schatz  [10],  or  Dryja  and  Widlund  [25]  for  a  proof  using  the  additive  Schwarz 
framework. 

We  remark  that  the  BPS  algorithm  and  the  algorithms  developed  in  this  chapter  are  it- 
erative substructuring  methods.  We  do  not  solve  subproblems  on  overlapping  subdomains 
(lij.  Instead,  we  solve  subproblems  on  non-overlapping  subdomains  fi,,  and  subproblems 
related  to  the  interface  r,j.  Estimates  for  the  convergence  rate  of  iterative  substructuring 
algorithms  are  independent  of  the  jumps  in  the  coefficients  of  the  elliptic  problems. 
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6.2      Iterative  Substructuring  Methods  for  the  Biharmonic 
Problem 

We  present  the  algorithms  and  the  analysis  only  for  the  bicubic  elements.  The  algorithms 
and  analysis  for  the  Argyris  and  Bell  elements  are  similar.  We  first  look  at  a  direct 
generalization  of  the  BPS  algorithm  and  show  that  the  condition  number  grows  at  least 
like  0{{H lhY)\  this  is  also  verified  by  numerical  experiments.  Modifications  are  needed 
to  get  well  conditioned  operators. 

As  in  the  additive  Schwarz  algorithms  cases,  we  use  Vq  and  Vq  to  denote  the  bicubic 
elements  associated  with  the  fine  and  coarse  triangulations,  respectively.  Vq  is  also  defined 
as  in  the  previous  chapter.  As  in  the  second  order  case,  we  need  the  subspaces  Vij,  defined 

by 

We  describe  several  possible  decompositions  of  the  finite  element  space.  Each  of  them 
results  in  an  algorithm.  Our  first  algorithm  resembles  one  of  the  algorithms  of  Bramble, 
Pasciak  and  Schatz  for  second  order  elliptic  problems. 

The  finite  element  space  Vq  can  be  represented  as: 

yQ  =  yQ^Yl^i:-  (61) 

Using  this  space  decomposition  and  corresponding  projections,  we  have 
Algorithm  6.2.1  Find  uy,  G  Vq  such  that 

with  an  appropriate  right  hand  side  g^ . 

We  recall  that  the  condition  number  of  the  BPS  algorithm  grows  like  (1  +  log^{H/h)). 
However,  for  the  biharmonic  operator,  the  condition  number  of  P^^^  grows  as  fast  as 
(1  +  \og{n/h)){n/h)^;  cf.  table  6.1  and  the  algorithm  is  not  very  practical.  If  we  use 
the  decomposition  6.1,  we  obtain  u/j  =  u//  +  J2ij  ^ij-  ^^  ^^^^e  that  the  only  choice  for 
ujf  is  the  interpolant  JlyHUh-  Because  of  the  appearance  of  the  second  derivatives  in  the 
interpolation  operator  UyH,  we  get  a  poor  bound  when  we  try  to  bound  ITlyHUhln''  by 
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|u/,|//3  and  this  is  reflected  in  the  poor  bound  in  the  condition  number  of  P^^h  This  is  a 
new  phenomenon  for  Hermitian  type  elements. 

One  way  to  overcome  this  difliculty  is  to  add  certain  vertex  subspaces  to  the  space 
decomposition  6.1.  We  then  have  more  freedom  to  choose  uh  and  can  avoid  the  inter- 
polation operator  UyH.    For  each  grid  point  Xk  6  A''(n),  we  define  a  vertex  space  V^, 

by 

Vj.  =  span{<;!l>f^}, 

where  <f>'^  is  the  basis  function  of  Vq  associated  with  the  mixed  second  derivative.   We 
then  have  a  space  decomposition 

Using  the  corresponding  projections,  we  obtain 
Algorithm  6.2.2  Find  Uh  £  V^  suck  that 

with  an  appropriate  right  hand  side  g^. 

In  each  iteration  of  this  algorithm,  we  have  three  sets  of  subproblems. 

1.  Compute  PyHVh  by  solving  Kyux  —  b.  This  is  a  finite  element  problem  in  Vq  . 
Here  KyH  is  the  stiffness  matrix,  and  dim(Vg^)  =  4A^,  where  N  is  the  number  of 
interior  coarse  grid  points. 

2.  Compute  Pv.jVh  by  solving  Kv,jX  =  b  exactly,  or  by  using  a  preconditioner.  This  is 
a  finite  element  problem  in  the  subdomain  n,j.  Techniques  for  constructing  precon- 
ditioners  for  the  problem  on  the  union  of  two  subdomains  can  be  used. 

3.  Compute  /V,  f/i  for  each  coarse  grid  point.  There  are  N  =  |A^^(n)|  such  problems. 
For  each  n,  Py,  Vh  is  computed  by  solving  KI^^x  =  6,  where  x  and  b  are  scalars 
and  K^^  —  a{4>'i^ ,  (j)^,^ )  is  the  1  by  1  stiffness  matrix  corresponding  to  the  degree 
of  freedom  associated  the  mixed  second  derivative  at  the  vertex  Xk-  Solving  such 
problems  are  trivial. 

74 


We  want  to  further  reduce  the  cost  per  iteration,  without  decreasing  the  rate  of  con- 
vergence. Using  the  fact  that  the  degrees  of  freedom  associated  with  the  second  order 
derivatives  are  only  weakly  coupled  to  the  other  degrees  of  freedom,  we  can  modify  Algo- 
rithm 6.2.2  slightly,  to  get  a  computationally  more  efficient  algorithm.  Let 

V;,  =  span{<AMQ|  =  0,l}, 

where  <^^,  \a\  =  0, 1,  are  the  basis  functions  of  Vq  associated  with  the  interior  vertices  of 
n,j.  We  note  that  Vjj  is  a  subspace  of  Vj-j,  with  vanishing  mixed  second  derivatives  at  all 
the  fine  grid  points  Xk  G  A''(n).  Using  Vij  and  Vq  instead  of  Vij  and  Vq  ,  we  obtain  the 
following  space  decomposition, 

v^  =  vq"+j:v„+  y:  v^,. 

*:eA''(n) 
Using  the  sum  of  the  corresponding  projections 

fceA''(n) 
We  obtain 

Algorithm  6.2.3  Find  u/,  6  V^  such  that 

P^^Kh  =  (Py^  +  Yl  Pv„  +      E      ^x.  >^  =  9k 

xteA'-Cfi) 

with  an  appropriate  right  hand  side  g^. 

We  note  that,  in  the  modified  algorithm,  we  reduce  the  size  of  the  the  coase  problem  and 
the  size  of  the  subproblems  associated  with  ft,j.  This  reduces  the  amount  of  work  per 
iteration  considerably.  We  need  to  solve  the  following  subproblems  in  each  iteration.* 

1.  Compute  PyHVh  by  solving  K^hX  =  b.  This  is  a  finite  element  problem  in  the  space 
Vq  .  We  note  that  KyH  is  a  principal  minor  of  KyH  and  dim(VQ')  =  ^dim(VQ^). 

2.  Compute  Py  Vh  by  solving  KyhX  =  b.  This  is  a  finite  element  problem  in  the 
subdomain  n,j  using  the  space  Vij.  This  task  is  similar  to  computing  Pv,jVh-  We 
note  however  that  Ky    is  a  principal  minor  of  A'v.^  and  dim(V,j)  =  ^dim(Vij). 
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3.  Compute  /V,  Vh  by  solving  Kv^  x  =  b;  cf.  Algorithm  6.2.2.  There  are  n  such 
problems,  where  n  is  the  total  number  of  interior  fine  grid  points.  Each  subproblem 
itself  is  similar  to  the  those  of  Algorithm  6.2.2.  Although  we  have  more  subproblems 
for  this  algorithm,  each  of  them  is  easy  to  solve.  The  total  amount  of  work  is 
proportional  to  the  number  of  unknowns. 

6.3      Condition  Number  Estimates 

We  have  the  following  theorem  for  the  operators  defined  above. 

Theorem  6.3.1    The  condition  numbers  of  the  operators  P^'\i  =  1,2,3,  for  the  additive 
algorithms  defined  above,  satisfy  the  following  estimates. 

Ac(p(*))   <  cin/h)\i-\-iog{n/h))  (6.2) 

K(P('))     <    Cil+log{H/h))\i  =  2,3.  (6.3) 

The  constant  C  is  independent  of  h  and  H . 

Before  we  prove  the  theorem,  let  us  discuss  the  decompositions  of  the  domain  and  the 
corresponding  partition  of  unity  for  the  iterative  substructuring  methods  constructed  in 
this  chapter. 

Let  x]j  and  xfj  be  the  two  end  points  of  r,j,  let  r  be  an  element  of  T^  and  let 

r(r)  =  min{dist(r,z]j),dist(r,3r?j)}. 

We  do  not  use  a  the  subscript  for  r{T),  since  we  always  work  on  only  one  substructure  fi,j 
at  a  time. 

First  we  consider  the  decomposition  • 

Associated  with  this  decomposition  of  the  domain,  there  is  a  partion  of  unity  {6,j}.  The 
properties  of  {0,j}  arc  summarized  in  the  following  lemma. 
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Lemma  6.3.1    There  exists  a  partition  of  unity  {Oij],  which  satisfies 

J29„{x)        =     I,  for  X  en,  (6.4) 

•J 

|^oIw".«>(t)     <     Cr{r)~*. 

If  we  use  additional  subdomains,  we  have  a  better  partition  of  unity  in  the  sense  that 
we  can  get  a  better  bound  for  the  W*'°°-norm  of  the  partition  functions.  We  consider  the 
following  decomposition  of  domain  il: 

Corresponding  to  this  decomposition  of  fi,  we  have  a  partition  of  unity  {0ij,Oic}-    The 
properties  of  this  partition  of  unity  are  summarized  in  the  following  lemma. 

Lemma  6.3.2    The  partition  of  unity  {^ij,^t}(i,i, /:  €  A'')  has  the  following  properties: 
LM^)+   E   ^t(x)      =       hforxed,  (6.5) 

ij  JteA" 

OijeC^iClij)     and    ^it  e  Co~(r(t)) 

|«.jk..oo(r)       <      C{min{r{T)-',h-']),     5  =  1,2, 
\^k\w-°°{T(xk))      ^      Ch~',     5  =  1,2. 
The  following  lemma  gives  bounds  of  the  spectrum  of  P^^h 
Lemma  6.3.3   The  operator  P^^^  satisfies  the  estimates: 

Ci{il  + login /h))in/h)'')-'a{uh,Uh)  <  aiP^'Kh,UH)  <  CMuh,UH),     Vu^  6  V^. 
The  constants  C\  and  C^  are  independent  of  h  and  H . 

Proof.  The  upper  bound  follows  from  the  finite  covering  property  of  {0,j}  and  prop- 
erties of  the  projections,  with  C2  <  A^c  +  1.  Here  iVc  =  4  for  rectangular  elements  and  3 
for  triangular  elements. 

To  prove  the  lower  bound,  we  need  to  get  a  partition  for  u/,  G  Vq.  Let  u//  =  JlyHU^ 
and  ly/i  =  u/i  —  u//.  Then 

Z?"t/;,(x,)  =  0,  |a|  =  0,l,     |^(^.)  =  0,     Vi,  6  A^(0).  (6.6) 
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Let  6ij  be  the  partition  of  unity  given  in  Lemma  6.3.1.  We  note  that  each  d,j  is  smooth 
except  at  the  two  substructure  vertices  common  to  Qi  and  Qj,  and  that  6,j  is  uniformly 
bounded.  Using  equation  (6.6),  we  can  show  that 

D''{e„w,){x,)  =  0,  \a\  =  0, 1,      ^-^f^i^,)  -  0,     Vx,  6  A^(n). 
It  is  not  difficult  to  see  that  the  interpolant  II v,(,(fl,jU>A)  is  well  defined  and  that  Uyh{0,jWh)  G 
Vij.  By  the  linearity  of  the  interpolation  operator  Uyh,  we  have, 

where  u//  G  V^  and  u,_,  G  V'.j.  We  now  estimate  |uw|//2  and  |u,j|//2.  By  Lemma  4.2.3,  we 
have 

From  Lemma  4.4.1,  we  obtain 

=    T1+T2  +  T3, 
where  6,(r)  =  |flij|tv..oo(i.).  Using  the  fact  that  \6,j\  <  1,  we  get 

Tcn.j 

We  have  the  estimates 

[H/h] 
E   /iVr(r)'  <  C  ^  1//  <  C(l  +  log(////i)),  (6.7) 


Thus, 


T2    =    C  Yl  0?(O|u'/.|'w.(r) 
Tcn.j 

Tcn., 

<    C(l+log(////i))|u.h|2vv„.„o(n.,). 
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We  note  that,  since  Wh  vanishes  at  the  substructure  vertices, 

Thus, 

Therefore, 

T,     =    C  Yl  Qlir)\wH\'mr)  <  C  E   hf^'^^f^lw^.-ia.,) 
Tcn.j  Tcii.j 

Thus, 

l«0lH^(n.,)  =  7^1  +  ^2  +  T3  <  C|Ti;;.|^,(n.^)  +  C(l  +  log(;7//i))|«;A|^^:,oo(n.,). 
From  Lemma  4.2.3,  we  know  that 

and 

Thus, 

l«ol^=(n,)  <  C((l  +  log(^//i))(^/ft)2|«,|2^,(jj_^j. 

Summing  over  all  fi.j  and  taking  the  finite  covering  property  of  Q.j  into  account,  we  get 

\^M\h(ti)  +  L  l«0l?f^(n.,)  <  C(l  +  log(^//i))(^//i)2|«,|2,,(j^j. 

The  lower  bound  now  follows  from  Lions'  lemma.  ■ 

Remark  6.3.1  For  any  decomposition  u^  =  Uff  +  ^.^  «,j,  the  coarse  function  Uff  has  to 
be  the  interpolant  UyHUh-  It  is  easy  to  construct  a  function  u/,  such  that 

Q 

Thus,  it  follows  from  Lemma  2.3.2  that  X„^n{P^^^)  <  CiU/h)-'^.  This  shows  that  the 
estimate  in  Lemma  6.3.3  is  almost  sharp.  This  bad  bound  for  the  condition  number  of 
P^^)  indicates  that  Algorithm  6.2.1  is  not  practical. 
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We  next  give  bounds  for  the  spectrum  of  f* '^^ 
Lemma  6.3.4    The  operator  P*^'  satisfies  the  estimates: 

Proof.  As  in  the  proof  of  the  previous  lemma,  the  upper  bound  follows  by  the  finite 
covering  property  of  {HijtTj^i,}  with  C2  =  Nc-¥  1. 

We  now  establish  the  lower  bound.   Let  «//  =  II^hU/,  6  Vq    and  u'/,  =  Uh  —  Uh-  By 
the  definition  of  JIw/zu/,,  we  know  that 

■(^*)  =  -nri'^k),      \o\  <  1     and     -K-^-{xk)  =  0. 


This  implies 


Let 


Then, 


dx"  dx""         '  dxdy 


^^(.,)  =  0,      H  <  1  and  ^x.)  =  ^x.). 


Ui^  =  nyh{6kWh)  and  u,j  =  UyhiOtjWh). 


"■''^'  =  ^'''  =  ^<''  =  S<"  =  »•  *'  ^  «"■' 


and 


«..(xO  =  ^(x.)-^(..)  =  o,    ^(..,  =  |!|„o,   v...A"(n, 


Therefore, 

-(ifc)<^r€  V^*,     VteA'^lfi)      and      u,jeV,j. 


^'*  "  ax%^ 


u,,  = 


By  the  linearity  of  Ili/h,  we  obtain 
0 


•i  *€A«(n) 
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We  now  estimate  |u,j|//2  and  |«n|//a.  As  in  the  proof  of  the  last  lemma,  we  have 


l«.ji//2(n.j)       =  ]C     l"'jl//2(r)   =      H     |nv^(^.j"'/.)|//2(r) 


<  c(5:0^(r)k,|^.(,)+  X]0?(r)|t.,|^,(,)+  ^0^(r)|u,,|i.(,A 

=    Ti  +  Tj  +  Ta 

<  C|u)A|^.(n,^)  +  C(l  +  log(/r//i))|u;;.|2^^„,oo(n.,) 

In  the  last  inequality,  we  have  used  the  following  inequalities  from  Lemma  4.2.4. 

and 

k'.l//^(n.,)  <  C(l  +log(^//»))l"/^l//^(n,)- 
By  Lemma  4.4.1  and  the  fact  that  |tfjt|w.,oo(T(xj^))  <  Ch~',  we  have 

l«rJ/fJ(n)       =       Wxk\]l^r(xu))  =  \^V^{^kWh)\]{2(r(x„)) 

<      C{  h-*\Wh\l7(r(^^))  +  /l"^k/.|//l(r(r,))  +  k/ilw^irCr,))  ^ 

In  the  last  step,  we  have  used  the  fact  that  D°Wh{xk)  =  0,  \a\  =  0, 1,  which  implies  that 

Wh\L^(r(x„))  <  Ch'^\wh\H^r{x,))     and     \wh\Hi{T{x,,))  <  CI^\^h\H^r(xk))- 

Combining  the  estimates  for  u,j  and  u^,  we  obtain 

which  by  Lions'  lemma  results  in  the  required  lower  bound  of  the  spectrum  of  P^^>.        I 
Lemma  6.3.5    The  operator  i'(^)  satisfies  the  estimates: 

Ci(l  +  log(^//i))-VuA,«0  <  a{P^^^Uh,Uh)  <  C2a{uh,uh),     Vu;.  e  V^. 
Proof.  The  proof  is  quite  similar  to  that  of  Lemma  6.3.4  I 
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6.4      Numerical  Results 

In  this  section,  we  report  on  some  numerical  results  with  the  iterative  substructuring 
methods  for  the  bihcirmonic  equation.  We  only  consider  the  bicubic  element  and  domain 
Q.  =  [0, 1]^.  As  for  the  additive  Schwarz  cases,  we  divide  the  fi  into  .V//  by  Nh  square 
subdomains  Qij.  Each  n,j  is  further  divided  into  elements.  Let  n\  be  the  total  number 
of  elements.  Then,  the  total  number  of  degrees  of  freedom  is  A(nh  -  1)^.  Let  H  =  j^ 
and  /i  =  -L.  We  stop  the  iteration  when  the  norm  of  the  residual  is  reduced  by  a  factor 
off  =  lO-''. 

In  the  first  set  of  experiments,  we  consider  the  iterative  substructuring  methods  with- 
out vertex  spaces;  cf.  algorithm  6.2.1.  To  study  the  dependence  of  the  condition  number 
of  F  on  (H/h)'^,  the  number  of  elements  inside  a  substructure,  we  restrict  our  experiments 
to  the  two  by  two  substructure  case.  From  table  6.1,  we  can  clearly  see  that  the  mini- 
mum eigenvalue  of  P  decreases  at  least  quadratically  with  H /h  and  that  the  maximum 
eigenvalue  of  P  remains  almost  fixed,  as  predicted  by  our  theory.  We  note,  however,  that 
although  the  condition  number  of  P  grows  very  fast,  the  number  of  iteration  grows  slower 
than  expected.  The  reason  might  be  that  convergence  rate  of  the  conjugate  gradient 
methods  depends  not  only  on  the  condition  number,  but  also  on  the  distribution  of  the 
eigenvalues  of  the  iteration  operator.  In  these  special  two  by  two  subdomain  cases,  we 
believe  that  only  a  few  eigenvalues  of  P  decrease  fast,  while  the  rest  of  the  eigenvalues 
remain  in  a  fixed  interval.  For  a  finite  difference  discretization,  a  similar  phenomenon  has 
been  observed  by  Chan,  E  and  Sun  [18].  They  also  gives  a  rigorous  proof.  We  believe 
however  that  this  is  not  true  for  problems  with  many  substructures. 

In  the  next  set  of  experiments,  we  consider  the  iterative  substructuring  methods,  with 
vertex  spaces,  i.e.  Algorithm  6.2.2.  For  a  comparison  with  last  set  of  experiments,  we 
also  restrict  our  experiments  to  two  by  two  substructure  cases.  From  table  6.2,  we  can 
see  that  the  maximum  eigenvalue  is  uniformly  bounded  by  three.  In  fact  the  maximum 
eigenvalue  decreases  slightly  as  H/h  increases.  The  minimum  eigenvalue  of  P  appears  to 
decrease  logarithmically  with  H/h,  the  number  of  elements  in  each  substructure. 
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total  #  unkns 

#  of  subdom 

-^min 

-^max 

K(F) 

#  of  iter. 

4x(4-l)^ 

1' 

0.99E-01 

2.1 

21 

9 

4x(8-l)^ 

2^ 

0.22E>01 

2.1 

95 

11 

4x(16-iy^ 

2^ 

0.53  E-02 

2.1 

400 

13 

4  X  (32 -1)2 

2^ 

0.13E-02 

2.1 

1600 

14 

Table  6.1:  ISM  Using  the  Bicubic  Element  V^^,  without  Vertex  Spaces 


total  #  unkns 

#  of  subdom 

'^min 

-^max 

k(P) 

#  of  iter. 

4x(4-iy^ 

2^ 

0.374 

2.1 

5.7 

8 

4x(8-iy2 

22 

0.222 

2.1 

9.6 

8 

4  X  (16 -1)2 

2^ 

0.146 

2.1 

14.5 

8 

4  X  (32 -1)2 

2^ 

0.103 

2.1 

20.5 

9 

4  X  (64 -1)2 

2^ 

0.0766 

2.1 

27.6 

9 

4  X  (100 -1)2 

2^ 

0.0645 

2.1 

32.6 

9 

Table  6.2:  ISM  Using  the  Bicubic  Element,  with  Vertex  Spaces 
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